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CHAPTER  1:  INTRODUCTION 


The  purpose  of  this  thesis  is  to  study  the  stability  characteristics  of  a  pair  of  coupled  van 
der  Pol  equations.  In  order  to  fulfill  this  purpose,  several  different  concepts  must  be 
discussed  prior  to  actually  looking  at  the  behavior  of  the  coupled  equations. 

The  first  issue  to  be  dealt  with  is  the  van  der  Pol  equation  itself;  this  topic  will  be  covered 
in  this  chapter  with  details  of  the  solution  to  the  van  der  Pol  equation  via  perturbation 
theory  given  in  Chapter  2. 

Chapter  3  will  explore  the  derivation,  existence  and  computation  of  Pade  approximants 
which  will  be  used  to  model  behaviors  of  the  van  der  Pol  equation  and  the  coupled 
equations.  If  the  reader  is  already  comfortable  with  the  characteristics  of  Pade 
approximants,  Chapter  3  can  be  skipped. 

In  Chapter  4,  the  model  for  the  coupled  van  der  Pol  oscillators  will  be  examined.  The 
scope  of  the  problem  of  characterizing  the  stability  behaviors  of  the  coupled  equations  will 
be  examined.  To  that  end,  the  Zero  Mean  Damping  (ZMD)  surface  will  be  introduced  and 
the  variational  equations  to  the  coupled  oscillators  will  be  derived.  Some  implications  of 
Floquet  theory  on  the  variational  equations  will  also  be  explored. 

Chapter  5  will  detail  the  solution  process  of  the  variational  equation  derived  in  Chapter  4. 
The  stability  transition  curves  along  the  ZMD  surface  will  be  introduced  and  the  method 
for  determining  them  will  be  presented. 

Chapter  6  will  summarize  the  results  of  Chapter  5  and  discuss  some  ongoing  work  and 
computational  difficulties  encountered  in  determining  the  stability  transition  curves. 
Further  studies  will  also  be  discussed. 


THE  VAN  DER  POL  EQUATION 
The  van  der  Pol  equation: 

u-e(l-u2)u  +  u  =  0  (1-1) 


was  introduced  by  Bismarck  van  der  Pol  in  a  paper  published  in  1927  [29,30]  to  model  the 
non-linear  resistance  of  certain  circuits  including  vacuum  tubes.  Since  that  time  Equation 
(1-1)  has  become  a  popular  mathematical  model  for  limit  cycle  oscillators. 

A  limit  cycle  oscillator,  is  an  oscillator  which  exhibits  self-sustained  periodic  behavior  It 
is  non-conservative,  due  to  a  damping  term,  and  non-linear  [2,7].  In  the  case  of  the  van 
der  Pol  equation,  the  non-linearity  is  also  due  to  the  damping  term.  Limit  cycle  behavior  is 
present  in  many  mechanical,  electrical,  and  biological  systems  [3,4,5,6].  Stoker[21]  and 
LaSalle[15]  showed  the  existence  of  a  unique  limit  cycle  for  Equation  (1-1)  for  all  e>0. 
The  limit  cycle  in  the  phase  plane  {u,  u'}  is  shown  in  Figure  (1-1)  for  the  case  of  6=0.9. 

The  shape  of  the  limit  cycle  for  the  van  der  Pol  equation  varies  markedly  with  the  value  of 
the  non-linear  parameter,  e.  For  e«1  the  limit  cycle  is  well  approximated  by  u=2  Cos(t), 
while  for  large  8  the  limit  cycle  remains  periodic  but  varies  significantly  from  sinusoidal 
[10,28,29].  For  large  8,  the  van  der  Pol  equation  can  be  used  to  model  relaxation 
oscillators  such  as  a  heart  beating,  a  clock  ticking  or  a  fish  swimming  [6,1 1,22,28]. 
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Figure  (1-1)  Phase-plane  for  VDP 
Oscillator  (e=0.9). 

Applications  of  relaxation  oscillators  is  the  motivation  for  this  thesis.  If  the  behavior  of  a 
chain  of  linked  van  der  Pol  equations  can  be  characterized,  it  may  be  possible  to  build  an 
electro-mechanical  system  which  uses  coupled  oscillators  to  produce  control  signals  for 
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locomotion.  As  an  example,  consider  a  pair  van  der  Pol  oscillators  in  the  relaxation  limit. 
Since  u(t)  is  characterized  by  a  rapid  transition  between  two  relatively  flat  states,  using  the 
signal  to  drive  a  pair  of  servomechanisms  may  be  a  very  inexpensive  way  to  drive  a  bipedal 
robot.  By  coupling  the  oscillators,  it  may  be  possible  to  produce  a  smooth,  responsive 
gait  with  little  or  no  computer  processing. 

Another  possible  application  for  coupled  van  der  Pol  oscillators  could  be  modeling  the 
synapses  in  a  fish's  spine  while  swimming  [6].  Here  the  chain  of  oscillators  may  be  30  or 
even  50  units  long  and  modeling  may  be  extremely  complicated. 

In  the  first  step  to  understanding  the  behavior  of  such  models,  it  is  desirable  to  find  the 
stability  characteristics  of  the  coupled  equations.  In  this  thesis,  the  scope  will  be  limited  to 
only  two  oscillators  which  could  be  modeled  as: 

x  -  e(1  -  x2  )x  +  x  =  eA(y  -  x)  +  eB(y  -  x) 

(1-2) 
y  -  s(l  -  y 2  )y  +  y  =  sA(x  -  y)  +  eB(x  -  y) 

where  sA  represents  the  displacement  coupling  and  eB  represents  the  velocity  coupling 
between  the  two  oscillators.  This  model  was  proposed  by  Rand  and  Holmes  [18]  and 
further  studied  by  Storti  and  Reinhall[  19,20,22,27].  In  future,  examining  the  stability  of  a 
chain  of  more  than  two  oscillators  could  be  done  by  extending  the  techniques  of  Storti  and 
Reinhall  and  this  thesis. 


CHARACTERIZING  THE  VAN  DER  POL  EQUATION 

The  first  step  in  understanding  the  stability  of  Equations  (1-2)  is  to  develop  the  solution  to 
the  van  der  Pol  equation.  Equation  (1-1)  has  been  studied  extensively  [1,  2,  7,  8]  and  a 
simple  approach  to  finding  solutions  using  Lindstedt's  method  [14]  is  presented  in  Chapter 

2. 


CHAPTER  2:  SOLUTION  TO  THE  VAN  DER  POL  EQUATION 


Perturbation  theory  is  one  of  the  most  powerful  techniques  used  in  solving  non-linear 
differential  equations  in  vibration  theory  [14,17].  In  concept,  perturbation  theory  is  very 
simple  and  can  be  explored  by  examining  the  van  der  Pol  equation  in  time  x: 

u-e(1-u2)u  +  u  =  0,        u(0)  =  0  (2-1). 

Using  Lindstedt's  method,  rescale  time  and  let  x  =  co  t.  Changing  variables,  Equation  (2- 
1)  becomes: 

u"  -  8(1  -  u2  )©  u'  +o) 2  u  =  0,        u'(0)  =  0  (2-2) 

where  the  primes  represent  differentiation  with  respect  to  the  rescaled  time  t.  Now  apply 
perturbation  theory,  by  considering  Equation  (2-2)  as  a  small  perturbation  of  the  linear 
system. 

u"  +  u  =  0,        u'(0)  =  0  (2-3) 

whose  behavior  is  well  understood.  It  is  reasonable  to  assume  that  u(t)  and  ©  may  be 
represented  by  power  series  in  s: 

u(t)  =  u0(t)  +  su1(t)  +  e2u2(t)  +  - 

G)  =  1  +  EC0,  +820)2  +•••  (2-4) 

Substituting  Equations  (2-4)  into  Equation  (2-2)  yields: 

Uq  +SU,"+82U2'  +  ---+(l+8(0j  +82C02  +---)2(U0  +E  U,  +E2U2  +•••)  = 

s{l-(u0  +  su,  +e2u2  +---)2}(l  +  e©,  +e2©2  +---)(Uo  +e  u|  +e2u2  +•••) 

(2-5) 

Collecting  like  powers  of  s,  the  first  few  equations  become: 


8°:      Uo+u0=0 

b1:      uj'+u,  =u;(l-u02)-2©,u0 

e2:      u2'  +  u2  =o1u;(l-u02)  +  u;(l-u02)-2u1u0u;-2co1u1  -2co2u0 -co,2u0 

(2-6). 

These  equations  are  then  solved,  in  order,  where  the  right  hand  side  is  completely  fixed  by 
the  solutions  to  the  previous  equations  and  determining  the  Oj's  and  any  undetermined 
coefficients  in  the  earlier  Uj(t)'s  to  suppress  secular  terms. 

Secular  terms  are  terms  that  if  not  canceled  in  some  manner,  would  give  solutions  to  the 
differential  equation  that  grow  with  time.  As  an  example,  consider  the  linear  Equation  (2- 
4)  with  an  added  driving  force: 

u"  +  u  =  Cos(t)  (2-7). 

Since  the  driving  force  matches  the  natural  frequency,  u(t)  will  grow  over  time  and 
become  unbounded;  the  solution  to  Equation  (2-7)  is  u(t)=l/2  t  sin(t),  which  clearly  grows 
over  time  In  applying  perturbation  theory,  the  goal  is  to  select  the  undetermined 
coefficients  in  Equation  (2-4)  to  eliminate  terms  which  match  the  frequency  of  Equation 
(2-3).  In  particular,  for  the  van  der  Pol  equation,  the  coefficients  of  any  Sin(t)  or  Cos(t) 
term  must  be  suppressed  when  they  appear  on  the  right  hand  sides  of  Equations  (2-6). 

As  a  demonstration,  consider  the  first  few  orders  of  e  in  Equations  (2-6).  Starting  with 
zeroeth  order,  the  solution  must  be  of  the  form: 

u0(t)  =  A0Cos(t)  +  B0Sin(t)  (2-8). 

Now  applying  the  initial  condition  in  Equation  (7),  B0=0  and  Ao  will  be  determined  while 
solving  the  first  order  equation.  Substituting  xo(t)  into  the  first  order  equation  and 
employing  trigonometric  identities: 

uj'  +  u,  =(A0Cos(t))'(l-(A0Cos(t))2)-2ca1(A0Cos(t)) 

(  \      \  1  (2-9) 

=  1 -A0  +  — AjJ  Sin(t)  +  -Aj  Sin(3t)-2©,A0Cos(t) 

In  order  to  suppress  the  secular  terms,  the  coefficients  of  Sin(t)  and  Cos(t)  must  be  set 
equal  to  zero,  there  are  three  possible  values  for  Ao  to  do  this,  Ao=-2,  0  or  2  and  (Oi  must 
be  zero.    Selecting  Ao=0  will  result  in  the  trivial  solution  for  the  zeroeth  order  equation, 
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but  either  2  or  -2  would  be  valid  and  give  solutions  to  Equation  (2-1)  with  opposite  signs. 
For  convenience,  take  the  solution  Ao=2.  The  right  hand  side  of  Equation  (2-9)is  now  2 
Sin(3t)  and  the  solution  is: 

u,(t)  =  B,Sin(t)--Sin(3t)  (2-10) 

and  imposing  the  initial  condition  of  Equation  (2-1)  again: 

u1(t)  =  |sin(t)-^Sin(3t).  (2-11) 

It  serves  little  point  to  continue  this  procedure  here,  but  it  is  worth  noting  that  both  the 
first  order  term  of  the  solution,  Ui(t),  and  the  first  order  term  in  the  frequency  expansion, 
coi,  were  determined  from  the  first  order  equation  of  Equations  (2-6).  This  pattern  will 
hold  while  solving  Equations  (2-6)  at  each  order  of  8.  Following  this  procedure,  it  is 
possible  to  find  solutions  for  all  u;(t)and  coj  up  to  a  high  value  of  i;  such  a  solution  is 
considered  an  i~  order  solution  to  the  differential  equation  [14]. 

This  can  be  a  tedious  solution  procedure,  but  the  main  advantage  is  that  any  degree  of 
accuracy  can  be  achieved.  In  particular,  if  the  solution  is  needed  to  some  high  order  n, 
this  method  will  allow  determining  the  solution  of  Equation  (2-1)  through  order  n  with 
sufficient  time  and  effort.  Fortunately,  this  method  (Lindstedt's  method)  can  be  realized 
using  a  programmable,  symbolic  algebra  package.  The  results  outlined  here  were 
determined  using  a  program  written  in  Mathematica™  [9,  32];  the  code  to  generate  these 
results  is  found  in  Appendix  A.  This  code  takes  advantage  of  some  additional  features  of 
the  van  der  Pol  equation  that  will  be  discussed  in  the  next  section. 

SYMMETRY  IN  THE  VAN  DER  POL  EQUATION 

Examining  Equation  (2-6)  with  symmetry  of  the  solution  in  mind,  it  can  be  seen  that  the 
right  hand  sides  of  the  equations  alternate  in  symmetry.  In  particular,  the  right  hand  side 
of  the  zeroeth  order  equation  is  0  which  is  even,  therefore  uo(t)  must  be  even,  which  it  is. 
The  right  hand  side  of  the  first  order  equation,  considering  the  fact  that  uo(t)  is  even,  must 
be  odd  resulting  in  an  odd  solution  for  ui(t).  This  pattern  continues  and  allows  the  rapid 
solving  of  each  order  of  the  power  series  solution. 

One  additional  property  of  the  equation  may  be  used  to  accelerate  the  solution  of 
Equations  (2-6).  Since  the  left  hand  side  of  every  order  in  8  is  the  same,  namely: 

LHS  =  u|'+u.  (2-12) 


then  the  homogeneous  solution  is  always  of  the  form: 

uu(t)  =  Ai  Cos(t)  +  B;  Sin(t)  (2-13) 

and  due  to  the  argument  of  the  above  paragraph,  Aj=0  if  i  is  odd  or  Bj=0  if  i  is  even. 

Even  more  fortunately,  the  particular  solution  for  each  order  in  e  is  also  predictable.  The 
right  hand  side  of  Equation(2-6),  after  suppressing  secularity,  will  always  take  the  form 
of: 


(2i+l),A=2 


RHS  =  a0  +     ^  ak  Cos(k  0  V  *  G  Even 

k=3 
(2i+l),A=2 

RHS=     £bkSin(kt)  VieOdd 


(2-14). 


k=3 


This  results  in  particular  solutions  of  the  form 

(2i+lXA=2 


xi,P(T)=     £     ,2k     Cos(kt)        VieEven 
k=3     k   - 1 

(2i+l),A=2       ^ 

x,p(t)=     X     rr^Sin^t)        VieOdd 


(2-15). 


The  overall  advantage  of  being  able  to  write  down  these  solutions,  is  that  Lindstedt's 
method  may  now  be  used  without  actually  solving  any  differential  equations.  This 
significantly  speeds  up  the  solution  generation  for  the  van  der  Pol  equation. 

After  removing  the  need  for  using  a  differential  equation  solver,  the  next  most  time- 
intensive  process  left  in  solving  each  order  of  the  van  der  Pol  equation  is  expanding  the 
right  hand  side  of  each  equation  and  using  the  trigonometric  identities  to  reduce  them  to 
the  form  of  Equation  (2-14).  One  very  simple  way  to  reduce  the  computational  effort  is  to 
use  Euler's  identity: 

ejt  =  Cos(t)  +  jSin(t),or 

_    ,  ,     ejt+e-jt  „  ,  ,     ejt+e-jt  (2-16) 

Cos(t)  = and    Sin(t)  = : 

2  2j 
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where  j  =  v^-1  Further,  in  using  Mathematica™,  the  expression  Exp(j  t)  is  very 
complicated,  so  it  is  advantageous  to  define  another  symbol,  say  Z  =  Exp(j  t),  and  then 
replace  each  sine  and  cosine  term  on  the  right  hand  side  of  Equations  (2-6)  with: 

/    v     Zk+Zk                 /    x     zk-zk 
Cos(kt)  = and     Sin(kt)  = (2-17) 

Using  these  expressions,  the  right  hand  side  equations  can  be  expanded  much  more  rapidly 
and  then  regrouped  back  into  terms  of  sines  and  cosines  to  implement  Equation  (2-15)  and 
suppress  secularity. 

The  first  few  terms  of  the  solution  of  Equation  (2-1)  are  presented  here: 

u  (t)  =  2  Cos  (t)  + 

/  3  Sin  (t)       1  v      2  /    Cos  (t)       3  5  N 

e Sin  (3  t)  +  e +  —  Cos  (3 1) Cos  (5 1)  + 

V      4  4  '        V        g  16  96  I 

,/     7Sin(t)        21  35  7  v  4 

M l77-+T77Sin(3t)-— -Sin(5t)+— -Sin(7t)  +0(e  ) 

v        256  256  576  576  /  (2-18) 

with  the  frequency  expansion  given  by: 

e2       lie4         35e6  678899e8  28160413e10  „ 

16       3072       884736       5096079360      2293235712000         v     ;  (2-19). 

Using  the  short-cuts  listed  above,  the  code  found  in  Appendix  A  was  developed.  The  first 
ten  terms  in  the  solution  of  the  van  der  Pol  equation  are  found  in  Appendix  C  along  with 
numerical  solutions  up  through  order  25.  The  solution  to  the  van  der  Pol  equation  was 
determined  analytically  (no  numerical  rounding)  up  through  order  75  using  this  code. 


CHAPTER  3:  BRIEF  DESCRIPTION  OF  PADE  APPROXIMANTS 


WHAT  IS  A  PADE  APPROXIMANT? 

A  Pade  approximant  is  a  rational  function  which  approximates  a  power  series  [12,13]  It 
is  represented  by  [L/P]  where  the  numerator  is  of  degree  L  and  the  denominator  is  of 
degree  P.  The  Pade  approximant  [L/P]  can  be  constructed  from  any  power  series  of 
degree  N: 

00 

f(z)  =  Zcizi  =c0+c1z  +  c2z2+---+cNzN+o(zN+1)  (3-1) 

i=0 

where  L  +  P  <  N,  and  takes  the  form: 

2  L 

r        -i     an  +  a,z  +  a-,z  +•••+&»  z  /  i^d^i\ 

[L / P]  =     °       '       . 2    2+        "  7v  +  0(zL+P+1 )  (3-2) 

b0+D,Z  +  D2Z    H hDpZ 

and  the  Pade  approximant' s  Maclaurin  series  expansion  must  match  Equation  (3-1) 
exactly  up  through  order  L  +  P.  In  order  to  make  [L/P]  unique  (since  multiplying  the 
numerator  and  denominator  by  any  value  would  create  another  approximant  with  the  same 
Maclaurin  series),  it  is  common  to  arbitrarily  set  bo=l  (or  divide  both  the  numerator  and 
denominator  by  bo),  thus  making  Equation  (3-2): 

2  L 

,        ,     an^"a1z  +  a■)z+••  '-("a,  z  /  , ,  D , .  \ 

tL/P1=  ut'    ,,'■,  +h  V  +Qz  <w> 

1  +  DJZ  +  D2Z    H hDpZ 

Note:  The  power  series  of  Equation  (3-1)  is  assumed  to  be  complex- 
valued  and  could  represent  a  Maclaurin  series  expansion  of  some  function 
f(z);  however,  the  power  series  need  not  be  complex-valued  in  order  to 
work  with  Pade  approximants  and  Equation  (3-1)  could  just  as  easily 
represent  the  Taylor's  series  expansion  of  a  real  valued  function,  or  an 
experimentally  obtained  polynomial  curve  fit,  or  a  power  series  solution 
obtained  with  perturbation  theory. 
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WHAT  ARE  PADE  APPROXIMANTS  FOR? 

There  are  three  main  reasons  for  constructing  Pade  approximants. 

The  first  reason  to  construct  [L/P]  is  to  produce  a  function  which  has  a  larger  radius  of 
convergence  than  the  original  power  series.  In  some  cases  [L/P]  may  have  a  significantly 
larger  radius  of  convergence  than  the  original  power  series,  or  possibly  even  converge  for 
all  values  of  z.  Even  if  the  approximant  does  not  have  a  larger  radius  of  convergence,  the 
approximant  may  yield  information  about  how  the  function  behaves  for  large  values  of  z. 
This  will  be  more  readily  seen  in  an  example. 

Example  3-1 

Consider  a  power  series: 

f(x)  =  l-x  +  x2-x3  +  --- 

which  is  the  Taylor's  series  for  1/(1  +  x)  has  a  radius  of  convergence  R<1. 
But  the  an  approximant  for  f(x)  is  given  by 


[1/1]  = 


1 


1  +  x 


which  is  the  original  function  and  is  finite  for  all  x  except  x=-land 
converges  for  R<oo.  One  can  see  the  behavior  of  f(x)  and  [1/1]  easily  in  a 
graph  (f(x)  is  the  dotted  line,  [1/1]  is  solid)  where  f(x)  includes  the  first  8 
terms  in  the  series: 


f(x)     &     [1/1] 


0 


Taylor's  Series 


[1/1]  Approximant 


12  3 

Figure  (3-1)  Taylor's  series  and  Pade  approximant  for  l/(l+x). 
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Here,  the  approximant  obviously  does  a  better  job  of  estimating  the 
function  for  x>l,  since  the  approximate  is  exactly  the  same  as  the  original 
function. 


The  second  reason  for  using  Pade  approximants  is  to  speed  up  the  evaluation  of  a  given 
power  series  when  working  inside  the  radius  of  convergence.  For  the  preceding  example, 
it  is  less  computationally  expensive  to  find  the  value  of  1+  x  and  its  reciprocal  than  it  is  to 
evaluate  the  Taylor's  series.  In  a  more  general  sense,  comparing  the  evaluation  of  two 
polynomials  of  degree  N/2  and  taking  their  ratio  is  a  less  costly  evaluation  than  calculating 
the  value  of  an  N  degree  polynomial,  due  to  the  expense  of  calculating  higher  powers  in 
the  N  degree  polynomial. 

The  third  reason  for  employing  Pade  approximants,  which  is  closely  related  to  the  first,  is 
to  more  accurately  model  the  asymptotic  behavior  of  a  function.  One  of  the  worst 
features  about  a  power  series  solution  is  that  the  highest  order  term  will  always  dominate 
when  x  is  large,  so  the  series  behaves  as  xN.  This  may  be  unfortunate,  especially  if  there  is 
some  indication  of  how  the  actual  solution  should  behave.  For  example,  if  the  function 
being  modeled  is  known  to  decay  as  1/x  or  is  known  to  approach  a  constant  value,  then  a 
Pade  approximant  should  be  more  successful  in  predicting  the  asymptotic  behavior  by 
careful  selection  of  L  and  P  when  constructing  the  [L/P]  approximant.  In  the  case  of  a 
function  which  should  decay  as  1/x,  selecting  P=L+1  will  probably  result  in  an 
approximant  that  also  behaves  as  1/x  in  the  limit  as  x  becomes  large  since  the  numerator 
will  behave  as  xL  and  the  denominator  as  xL+1.  In  the  case  where  the  function  should 
approach  a  constant,  picking  L=P  will  most  likely  produce  an  approximant  with  the 
desired  behavior.  (Example  3  below  will  explore  some  of  this  more  clearly.) 

CALCULATING  PADE  APPROXIMANTS 


Procedure 

Calculating  a  Pade  approximant  is  a  straight-forward  process  [13].    Starting  by  setting 
Equation  (3-1)  equal  to  Equation  (3-3): 

j  n      a0  +  a,  z  +  a,  z2  +•  •  H-a,  zL 

c0+Clz  +  c2z2+.~+cNzN  =    °       '         2  L  (3-4) 

l  +  b,Z  +  D2Z    H hDpZ 

Next,  multiply  both  sides  by  the  denominator  of  the  right  side: 

(c0  +  c,  z  +  c2  z2  +•  •  -+cN  zN  )(l  +  b,  z  +  b2  z2  +•  •  -+bp  zp)  = 


a0  +  a,  z  +  a2  z2  +•  •  -+aL  zL 


(3-5) 
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Multiplying  Equation  (3-5)  out  and  collecting  like  powers  of  z,  yields: 


CL-P+1        CL-P+2  CL 

CL-P+2        CL-P+3        ""         CL+1 


CL  CL-P+1 


'L+P-l 


"b," 

cL+1' 

bp  , 

= — 

CL+2 

U . 

_CL+P_ 

(3-6) 


where  Cj  =  0  if  j<0  for  consistency,  and 


ao  -  co 


a,  -c,  +  b,c0 
a2  =c2  +b,c,  +b2c0 

min(L,P) 

aL=cL+    ZbicL-i 

i=l 


(3-7). 


The  linear  system  in  Equation  (3-6)  may  be  solved  for  the  b;'s  by  any  suitable  method,  and 
then  the  b;'s  used  in  Equations  (3-7)  to  find  the  ai's.  Following  this  procedure,  any  order 
Pade  approximant  can  be  constructed  for  L  +  P  <  N. 

Equations  (3-6)  and  (3-7)  can  be  easily  coded  into  any  programming  language  or  written 
as  functions  in  many  mathematics  packages. 

EXAMPLES 


Example  3-2 

Find  [4/4]  for  the  Taylor's  series  approximation  to  Sin(x).    Starting  with 
the  Taylor's  series: 


3  5  7 

XXX 


^T+o(x>) 


and  [4/4]  is 


[4/4] 


a0  +a,x  +  a2x2  +a3x3  +a4x4 
l  +  b,x  +  b2x2  +b3x3  +b4x4 


Using  Equation  (3-6): 
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1             0  -1/3!  0 

0  -1/3!        0  1/5! 

1/3!        0          1/5!  0 

0          1/5!          0  -1/7! 


p>4] 

"  1/5!  " 

p>4] 

"11/5880" 

b, 

0 

b, 

0 

b2 

— 

-1/7! 

=> 

b2 

— 

3/49 

Lb,_ 

0 

Lb,_ 

0 

a0=0 
a,  =1  +  0  =  1 

a2  =0  +  0  +  0  =  0 
a3  =-l/6  +  0  +  3/49  +  0  =  -31/294 


0+0+0+0+0=0 


so  the  final  solution  is: 


x  - 


[4/4]  = 


31 
294 


3_  2        11 
1  +  49X   +5880: 


+ 


O(x') 


Note.  Pade  approximants  do  not  perform  much  better  than  regular  power 
series  in  predicting  periodic  functions,  this  example  was  included  for 
numerical  interest  as  will  be  discussed  later. 


Example  3-3 


Another  example  of  finding  a  Pade  approximant  demonstrates  how 
asymptotic  behavior  for  Pade  approximants  may  be  significantly  better  than 
a  normal  power  series    Consider  the  function: 


f(x)  = 


V? 


with  Pade  approximant  [3/4] 


2 


+  x 


\+xl+-x 
8 


Selecting  a  [3/4]  approximant  increases  the  likelihood  of  correctly 
modeling  the  1/x  behavior  of  f(x)  for  large  x.  Attempting  to  calculate  the 
[3/4]  approximant  yields  a  [2/4]  approximant  which  shows  that  not  all 
degree  approximants  are  necessarily  unique.  Below  is  a  graph  of  f(x)  (solid 
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line),  [3/4]  (dotted  line)  and  the  eighth  order  Taylor's  series  for  f(x)  (dash- 
dot).  In  this  case,  the  [3/4]  approximant  does  a  significantly  better  job  of 
modeling  the  original  equation  and  tends  to  zero  (as  1/x2)  for  large  x. 


-l 


6  8  10 

[3/4]  Approximant 


Taylor's  Series 


Figure  (3-2)  Taylor's  series  and  a  [3/4]  approximant  to  f(x) 


Vf 


+  X' 


STABILITY  OF  PADE  APPROXIMANTS 


Condition  Number  Considerations 

Following  the  procedure  set  out  here  for  finding  a  Pade  approximant  includes  solving  the 
linear  system  represented  by  Equation  (3-7)  and  Equations  (3-7).  While  Equations  (3-7) 
present  no  difficulties  computationally,  the  solution  to  Equation  (3-6)  may  be  subject  to 
significant  loss  of  numerical  accuracy.  This  is  of  particular  concern  when  combined  with 
the  fact  that  the  accuracy  of  Pade  approximants  is  highly  sensitive  to  the  differences 
between  coefficients  in  the  approximant.  (The  need  for  very  high  numerical  accuracy  is 
discussed  in  detail  in  Reference  [13],  section  2. 1 .) 

To  illustrate  how  significant  the  numerical  results  may  be  affected  when  solving  Equation 
(3-6),  one  can  examine  the  condition  number  (or  a  reasonable  approximation  for  the 
condition  number)  of  the  matrix  C: 


CL-P+1        CL-P+2 
CL-P+2        CL-P+3 


'L+P+l 
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'L+l 


'L+P-l 


(3-8). 


To  gain  some  measure  of  the  magnitude  of  the  condition  number  of  this  matrix,  consider 
the  calculations  carried  out  in  Example  3-2.  Below  is  a  table  which  shows  how  the 
condition  number  of  C  varies  with  the  degree  of  the  Pade  approximant  used  to  estimate 
the  function  f(x)=Sin(x)  and  the  number  of  lost  digits  in  accuracy  is  assumed  to  be 
Logio(condition  number): 


L=P      Cond  # 


Lost  Digits 


4 

4.72  104 

4 

5 

1.66  105 

5 

6 

4.27  107 

8 

7 

3.40  109 

10 

8 

1.74  1012 

12 

9 

2.46  1014 

14 

10 

2.08  1017 
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Table  (3-1)  Condition  number  and  loss  of 
accuracy  for  Example  (3-2). 


The  importance  of  finding  an  accurate  method  to  solve  Equation  (3-6)  is  apparent  with 
respect  to  these  results,  and  it  would  be  wise  to  check  the  condition  number  of  C  when 
employing  Pade  approximants. 

There  is  one  obvious  way  to  avoid  the  numerical  problems  discussed  here.  By  employing 
one  of  the  powerful  math  packages  available  which  can  do  exact  arithmetic  such  as 
Mathematica™  or  Maple™,  Equation  (3-6)  can  be  solved  with  no  loss  of  accuracy  when 
the  coefficients  Cj  are  known  as  precise  fractions.  While  it  may  seem  optimistic  to  want  to 
know  the  Cj's  exactly,  in  many  applications  where  Pade  approximants  are  used,  the  power 
series  being  replaced  is  generated  exactly  from  a  technique  such  as  perturbation  theory.  In 
particular,  in  non-linear  vibration  analysis,  the  solution  to  a  differential  equation  may  be  a 
power  series  in  the  non-linear  parameter,  8,  known  to  50  or  even  100  terms  which  are 
calculated  with  fractions  and  known  precisely.  This  will  be  discussed  more  fully  when 
examining  applications  of  Pade  approximants  to  vibration  analysis. 
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What  is  a  defect? 

In  simple  terms,  a  defect  is  a  failure  of  a  Pade  approximant  to  accurately  approximate  a 
power  series  in  the  neighborhood  of  a  root  of  the  denominator  of  the  approximant.  Or  in 
mathematical  terms,  a  defect  occurs  around  x=xD  when  l+b,xD +b2Xo+ •• +bPXD  =0. 
Borrowing  on  the  language  of  complex  analysis,  the  defect  occurs  where  [L/P]  has  a  pole 
at  xd  or  zD. 

How  to  handle  a  defect 

The  presence  of  the  pole  in  [L/P]  may  or  may  not  be  a  result  of  numerical  inaccuracy  in 
calculating  the  approximant,  and  on  a  first  look  it  seems  reasonable  that  there  would  be 
cases  where  we  would  expect  the  poles  to  be  present.  The  most  obvious  reason  for  the 
[L/P]  approximant  to  have  a  pole  would  be  if  the  original  function,  f(z),  also  has  a  pole. 
In  this  case,  the  behavior  of  the  approximant  will  depend  on  the  nature  of  the  pole  in  f(z), 
if  the  pole  is  an  essential  pole  (there  is  no  finite  n  such  that  (z-zD)n  f(z)  is  analytic  at  zD) 
then  the  approximant  will  always  have  a  defect  at  zD,  regardless  of  the  degree  of  [L/P] 
However,  if  the  pole  in  f(z)  is  pole  of  finite  order,  it  may  be  possible  to  remove  the  defect 
in  the  approximant  by  constructing  [L/P]  with  L  and  P  large  enough.  There  is  no 
guarantee  that  such  an  approximant  is  possible  to  construct,  but  taking  higher  degree 
approximants  frequently  removes  defects  when  applied  to  perturbation  theory. 

A  defect  may  also  occur  in  an  approximant  when  f(z)  has  no  poles.  In  this  case,  the  defect 
could  occur  due  to  computational  errors  which  can  be  significant  as  discussed  in  the  last 
section.  To  handle  defects  of  this  nature,  there  are  two  basic  approaches.  One  approach 
is  to  use  a  more  accurate  method  to  compute  the  solution  of  Equation  (3-6).  This  may 
entail  finding  the  coefficients  in  the  original  power  series  to  greater  accuracy,  using  a 
machine  capable  of  higher  accuracy  or  resorting  to  rational  approximants  of  the 
coefficients  of  the  original  power  series  and  then  solving  Equation  (3-6)  exactly.  The 
second  method  would  be  to  try  calculating  the  approximant  to  lower  degree  in  the  hopes 
that  the  C  matrix  condition  number  will  be  enough  smaller  to  prevent  the  resulting 
approximant  from  having  a  defect. 

Another  cause  for  a  defect  to  occur  is  trying  to  approximate  a  complicated  function  by  too 
low  an  order  Pade  approximant.  If  this  is  the  case,  merely  finding  a  higher  order 
approximant  will  remove  the  defect  with  little  effort.  Coupling  this  comment  with  the  last 
point  in  the  previous  paragraph,  it  seems  that  trying  various  degree  approximants  (higher 
and  lower)  until  finding  one  that  is  defect-free  is  a  reasonable  approach,  however, 
checking  the  condition  number  of  C  is  an  intelligent  starting  point  to  determine  if  the 
difficulty  is  computational. 

More  information  on  making  the  calculation  of  Pade  approximants  numerically  stable  can 
be  found  in  Reference  [13]. 
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APPLYING  PADE  APPROXIMANTS  TO  PERTURBATION  RESULTS 

Having  briefly  explored  one  type  of  perturbation  theory  application,  there  are  some 
obvious  opportunities  to  apply  Pade  approximants  to  good  ends. 

One  of  the  most  direct  applications  for  Pade  approximants  in  vibrations  is  to  predict  the 
frequency  response  of  a  system  in  terms  of  the  non-linearity  parameter  8.  There  are  many 
different  methods  for  estimating  the  natural  frequency  of  slightly  non-linear  systems  (8  « 
1),  but  many  are  inaccurate  if  the  system  is  tuned  such  that  s  is  slightly  larger.  Given  that, 
it  is  reasonable  to  calculate  a  Pade  approximant  for  the  frequency  of  a  non-linear  system 
obtained  with  a  perturbation  technique  in  an  effort  to  accurately  predict  the  response 
frequency  for  larger  s.  As  a  starting  point,  it  is  assured  that  a  straight  power  series  for  co 
will  fail  to  accurately  predict  the  frequency  an  oscillatory  solution  as  s  increases  except  in 
very  special  circumstances.  In  the  general  case  the  frequency  will  not  grow  as  en  which  is 
the  asymptotic  behavior  of  the  power  series  for  large  s.  This  can  be  seen  by  examining  the 
power  series  expansion  for  the  frequency  of  the  van  der  Pol  equation. 

Example  3-4 

Returning  to  the  van  der  Pol  oscillator,  Equation  (1-1),  and  using  Lindstedt's  method,  the 
frequency  can  be  found  to  high  order  in  e.  (This  result  was  derived  in  Chapter  2.)  The 
first  few  terms  are: 

,      1     2        17      .  35        ,         678899 

co  =  1-— e2  + e4  + e6 s8  +0(e10) 

16         3072  884736         5096079360  v      ' 

Taking  the  [4/4]  approximant: 

492389  ,   433361 

1  + e2+ s4 

r4  /  4i  =   2698560    34541568 
1   J     661049  ,   1285087 

1  + 8    + 8 

2698560         57569280 
Graphing  both  the  power  series  for  co  (dashed)  and  [4/4]  (solid): 
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CO 


Figure  (3-3)  Power  series  and  Pade  approximant  to 
the  frequency  of  a  van  der  Pol  limit  cycle. 


From  experimental  data,  the  van  der  Pol  oscillator  is  in  fact  a  softening  system,  meaning 
the  frequency  of  the  response  decreases  with  increased  non-linear  effects.  Also,  just  from 
examining  the  graph,  the  behavior  of  the  [4/4]  approximant  seems  to  continue  the  pattern 
of  the  power  series  expansion,  before  8  became  large  enough  to  make  co  start  to  grow 
quickly.  Both  point  to  the  fact  that  the  Pade  approximant  does  a  good  job  of  extending 
the  approximation  for  frequency  response  as  8  increases. 

Large  s  Concerns 

Another  area  where  Pade  approximants  are  very  useful  in  vibration  theory  is  in  application 
to  relaxation  oscillators.  Relaxation  oscillators  have  numerous  applications  in  biological 
systems  including  how  the  heart  beats,  how  fish  swim  (modeling  the  synapses  in  the  fish's 
central  nervous  system),  and  have  possible  application  as  control  circuits  for  robotics.  The 
relaxation  oscillator  is  modeled  by  a  non-linear  differential  equation  where  the  non- 
linearity  is  allowed  to  become  large  compared  to  the  linear  components  in  the  system, 
meaning  in  the  neighborhood  of  10  or  100  times  larger.  Consider  the  relaxation  limit  for 
the  van  der  Pol  oscillator,  Equation  (1-1),  which  behaves  as  a  relaxation  for  large  8.  As 
seen  in  Example  3-4,  the  Pade  approximant  for  the  frequency  more  reasonably  predicts 
behavior  in  the  relaxation  limit. 


Application  to  Coupled  van  der  Pol  Oscillators 

The  main  reason  for  presenting  this  information  on  Pade  approximants  is  in  applying  the 
theory  to  the  stability  transition  curves  that  will  be  derived  in  Chapter  5.  For  the 
remainder  of  this  thesis,  some  basic  familiarity  with  Pade  approximants  will  be  assumed  as 
they  will  be  used  to  estimate  the  frequency  of  the  limit  cycle  of  Equation  (1-1),  the  Zero 
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Mean  Damping  Surface  described  in  Chapter  4,  and  the  Stability  Transition  Curves  found 
in  Chapter  5. 


CHAPTER  4:  COUPLED  OSCILLATOR  ANALYSIS 


PROPOSED  MODEL  FOR  COUPLED  OSCILLATORS 

As  discussed  in  the  introduction,  the  pair  of  van  der  Pol  oscillators  (x  and  y)  will  be 
modeled  by  two  van  der  Pol  equations  coupled  with  a  linear  spring  and  a  linear  damper. 
The  coupling  will  be  "weak"  in  the  respect  that  the  spring  and  damper  will  both  be  order 
s,  in  particular  the  stiffness  constant  will  be  8 A  and  the  damping  constant  will  be  eB  where 
8  is  the  same  non-linear  parameter  as  in  the  original  oscillator.  The  coupled  equations  take 
the  following  form  [18]: 

x  -  e(  1  -  x2  )x  +  x  =  eA(y  -  x)  +  eB(y  -  x) 

(4-1) 
y-e(l-y2)y  +  y  =  eA(x-y)  +  eB(x-y) 

where  dots  represent  derivatives  with  respect  to  time  x. 

The  object  is  to  determine  the  existence  and  stability  of  solutions  to  Equation  (4-1).  In 
particular,  it  is  desired  to  find  what  values  for  eA  and  eB  result  in  orbitally  stable 
solutions. 

In  this  context,  orbitally  stable  is  most  easily  described  in  the  phase  plane.  Consider  the 
phase  plane  of  the  original  van  der  Pol  equation  plotted  in  Figure  (1-1);  the  limit  cycle 
presented  is  orbitally  stable  In  mathematical  parlance,  for  any  small  perturbation  of  the 
system  away  from  the  described  limit  cycle,  there  exists  a  finite  distance  8  such  that  the 
maximum  distance  between  any  point  on  the  disturbed  path  and  the  original  limit  cycle  is 
less  than  or  equal  to  8  [14]  Understanding  the  concept  of  orbital  stability  is  essential  to 
comprehending  the  behavior  of  coupled  oscillators. 

By  inspection,  an  exact  in-phase  motion  of  the  coupled  oscillators  is  possible  with 
x(x)=y(x).  In  that  case,  the  coupling  terms  identically  vanish  and  x(x)  and  y(x)  will  have 
to  satisfy  the  original  van  der  Pol  equation  as  found  in  Chapter  1  with  x(x)  =  y(x)  =  u(x). 
In  this  situation,  the  behavior  of  each  of  the  coupled  oscillators  is  exactly  the  same  as  the 
original  solution,  namely  each  oscillator  would  exhibit  an  orbitally  stable  limit  cycle  in  its 
own  two  dimensional  phase  plane,  as  seen  in  Figure  (1-1). 
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In  the  same  vein,  it  is  possible  to  take  advantage  of  the  known  properties  of  u(x)  by 
considering  x(x)  and  y(x)  as  functions  which  vary  from  u(x)  by  small  amounts.  Treating  x 
and  y  as  variations  away  from  the  solution  of  the  van  der  Pol  equation  clarifies  the 
meaning  of  "stability"  as  applied  to  the  coupled  equations  of  Equation  (4-1).  While 
studying  the  variations  away  from  u(x),  determining  the  stability  of  x  and  y  is  equivalent  to 
studying  the  origin  in  the  variational  space.  This  will  be  addressed  in  more  detail  after 
deriving  the  variational  equations. 


DERIVATION  OF  THE  VARIATIONAL  EQUATIONS 

Treating  x(x)  and  y(x)  as  small  variations  from  u(x),  the  stability  of  the  solutions  of 
Equation  (4-1)  can  be  determined  by  characterizing  the  behavior  of  the  variations. 
Consider  both  functions  as  small  variations  from  the  limit  cycle  u(x)  [28]: 

£  =  x-u  =>  \  =  t  +  u 

(4-2) 
rj  =  y-u  =>  y  =  T]  +  u 

where  t,  is  the  variation  of  x(x)  from  u(x)  and  r\  is  the  variation  of  y(x)  from  u(x). 
Replacing  x  and  y  in  Equation  (4-1)  and  rearranging  terms  yields: 

£-s(i-£2-2u£-u2)£  +  (1  +  28uu)£  +  {u-e(1-u2)u  +  u}  = 

bA(ti-S)+eB(t}-4) 

r|  -  s(l  -  r)2  -  2ur]  -  u2  )t]  +  (l  +  2euu)r)  +  {ii  -  e(l  -  u2  )u  +  u}  = 

8A(4-ti)  +  8b(S-ti) 


(4-3) 


where  the  terms  in  braces  are  the  original  van  der  Pol  equation  and  are  identically  equal  to 
zero,  resulting  in: 

l  -  s(l  -  £2  -  2u£  -  u2)£  +  (1  +  2euu)£  =  bA(ti  -  $)  +  £b(ti  - i) 
r\  -  e(l  -  r)2  -  2uti  -  u2  )r)  +  (1  +  2suu)r|  =  sA(^  -  r\)  +  sByi  -  fi) 

In  the  interest  of  studying  only  small  variations  from  the  limit  cycle,  non-linear  terms  in  £ 
and  r|  are  neglected: 
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l  -  e(l  -  u2 )i  +  (1  +  2suii)$  =  sA(ti  -  I;)  +  eb(ti  -  $) 

r|  -  s(l  -  u2  )ti  +  (1  +  2euu)r|  =  eA(^  -  ti)  +  eb(£  -  rj) 

In  order  to  simplify  the  coupling  terms  appearing  on  the  right  hand  sides  of  Equation  (4- 
5),  two  new  equations  can  be  formed  by  adding  and  subtracting  the  two  Equations  (4-5) 
and  letting: 

w  =  £  +  T| 

•      -  (4-6). 

So  Equations  (4-5)  become: 

w  -  e(l  -  u2  )w  +  (l  +  2euu) w  =  0 

/  x  (4"7> 

v  -  e(1  -  2B  -  u2  )v  +  (1  +  2suu  +  2sA)v  =  0 

By  inspecting  Equations  (4-7)  it  can  be  seen  that  both  equations  are  exact  differentials 
which  can  be  integrated  with  respect  to  time  to  yield: 

w-sll-  u2)w  +  w  =  k, 

/  x  (4~8)- 

v  -  e(l  -  2B  -  u2  )v  +  (1  +  2eA)v  =  k2 

While  studying  stability,  the  solutions  for  v(i)  and  w(x)  can  be  recast  to  further  simplify 
Equation  (4-6).    The  particular  solutions  to  Equations  (4-8)  are  constants,  in  particular: 

wp  =k, 

k2  (4-9). 

Vp  ~  l  +  2eA 

These  solutions  do  not  vary  with  time  and  can  be  removed  by  a  simple  coordinate 
transform  in  the  variational  space:  w-»  w-wp  and  v-»v-vp.  Such  a  transformation  will  not 
alter  the  stability  of  the  solutions,  since  stability  is  controlled  by  the  exponential  behavior 
of  the  functions.  Then  the  varational  equations  become: 

w  -  e(l  -  u2  )w  +  w  =  0 

:  .  (4-10). 

v-e(l-2B-u2)v  +  (l  +  2eA)v  =  0 
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In  order  to  study  the  stability  of  the  coupled  van  der  Pol  oscillators,  both  of  Equations  (4- 
10)  must  be  explored  and  characterized.  In  order  to  discover  the  nature  of  the  first 
equation  in  Equation  (4-10),  consider  a  small  variation  from  a  single  van  der  Pol  solution, 
u(x): 

z(t)  =  u(t)  +  $(t)  (4-11) 

where  u(t)  is  a  solution  of  the  van  der  Pol  equation,  Equation  (1-1).    Substituting  z  into 
Equation  (1-1)  results  in: 

u  +  £  +  e(l-(u  +  £)2)(u  +  £)  +  u  +  £  =  0  (4-12) 

which  when  expanded  and  eliminating  the  higher  order  terms  in  £  results  in: 

{u  +  s(1-u2)u  +  u}+£+e(i-u2)£  +  £  =  0  (4-13) 

where  the  term  in  braces  is  again  the  van  der  Pol  equation  and  is  identically  zero,  resulting 


in: 


£+e(1-u2)£  +  £  =  0  (4-14). 

Equation  (4-14)  is  the  behavior  of  a  small  variation,  £  away  from  the  orbitally  stable  van 
der  Pol  limit  cycle,  therefore  £  must  be  orbitally  stable.  Comparing  Equation  (4-14)  to  the 
first  of  Equation  (4-10),  w  and  £  have  identical  solutions,  therefore  w(x)  must  be  orbitally 
stable. 

Understanding  the  stability  of  the  first  Equation  (4-10)  leaves  only  the  behavior  of  the 
second  variational  equation: 

v-e(l-2B-u2)v  +  (l  +  2eA)v  =  0  (4-15) 

to  be  characterized  in  order  to  understand  the  stability  of  Equations  (4-1).  Equation  (4- 
15)  will  be  referred  to  as  the  variational  equation  for  the  remainder  of  this  thesis. 

SIGNIFICANCE  OF  THE  ZERO  MEAN  DAMPING  SURFACE 

Studying  the  stability  of  the  variational  equation  requires  characterizing  the  solutions  of 
Equation  (4-15)  in  the  three  dimensional  parameter  space  {  e  ,  8 A,  sB  }.  Ultimately,  it 
should  be  possible  to  find  the  surface  in  {  e  ,  sA,  eB  }  space  which  corresponds  to  the 
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transition  from  stable  to  unstable  oscillations  of  the  coupled  equations.  Such  an  analysis  is 
feasible,  however  it  would  be  formidable  to  complete  using  analytical  methods.  Such  an 
analysis  using  numerical  techniques  is  underway  by  Reinhall  and  Storti[  19,20]. 

In  order  to  reduce  the  complexity  of  the  problem,  the  parameter  space  may  be  reduced  to 
two  dimensions  by  fixing  or  specifying  any  of  the  three  parameters.  Such  an  analysis  for 
the  case  of  B=0,  i.e.  no  damping  between  the  coupled  oscillators,  is  presented  by  Storti 
and  Reinhall  [26]. 

Another  obvious  two  dimensional  subspace  for  studying  the  variational  equation  is  the 
Zero  Mean  Damping  (ZMD)  surface  of  the  parameter  space.  The  ZMD  surface  is  the 
locus  of  values  for  the  damping  parameter  eB,  as  a  function  of  e,  resulting  in  zero  time 
average  damping  in  the  variational  equation.  The  ZMD  surface  is  a  logical  place  to  study 
the  behavior  Equation  (4-15)  because  it  marks  the  boundary  between  the  region  of 
negative  mean  damping,  where  all  solutions  to  Equation  (4-15)  will  grow  over  time  due  to 
the  energy  added  by  negative  damping,  from  the  region  where  damping  is  positive  on 
average  and  stable  solutions  may  be  possible.  By  examining  Equation  (4-9),  the  parameter 
sA  does  not  appear  in  the  damping  term,  so  the  ZMD  surface  will  be  a  function  of  e  and 
eB  only.  The  ZMD  surface  is  found  by  requiring  the  time  average  of  the  coefficient  of  v 
to  vanish  in  Equation  (4-15): 

(l-2BZMD-u2)-0  (4-16) 

where  the  angled  brackets  indicate  time  averaging.  Clearly  BZmd  will  be  a  function  of  e 
since  u(t),  found  in  Chapter  1,  is  a  power  series  in  e. 

The  balance  of  this  thesis  will  be  dedicated  to  characterizing  the  behavior  of  the  variational 
equation,  Equation  (4-15),  along  the  ZMD  surface. 

CALCULATING  THE  ZMD  SURFACE 

Performing  the  time  averaging  using  the  solution  for  u(t)  found  in  Chapter  1  and  is  a 
tedious,  but  straight-forward  process.  Since  BZMD  is  not  a  function  of  time,  Equation  (4- 
16)  can  be  rewritten  as: 

Bzmd=~(u2)  (4-17). 

Using  the  solution  u(t)  generated  in  Chapter  1 ,  the  ZMD  surface  can  be  calculated  to  the 
same  order  in  s  as  u(t).  The  first  few  terms  of  BZMd  are: 
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B 


ZMD 


1  11  lle    4859s6    12921629s8 

2  ~  32  +  9216  + 10616832  ~ 152882380800 


0(s'°)      (4-18). 


A  Pade  approximant  to  the  ZMD  surface  through  O(s20)is  given  by: 

-0.5  -  0.18562  -  0.0669e4  -  0.01 10e6  -  0.00132 e8  -  0. 0000589c 10 


BzMD  = 


1.  +  0.307  e2  +  0.1  He4  +  0.0163  e6  +  0.00202e8  +0. 0000654  e10      (4-19) 


where  the  fractions  which  actually  appear  in  the  Pade  approximant  have  been  replaced 
with  3  decimals  of  precision.  Figure  (4-1)  is  a  plot  of  the  Pade  approximant  [10/10]  of 
Bzmd  calculated  using  exact  coefficients. 
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Figure  (4-1)  Pade  approximant  to  the 
Bzmd  vs.  s. 
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Figure  (4-2)  Pade  approximant  to  the 
ZMD  surface  in  3-D  parameter  space. 


Figure  (4-2)  shows  no  additional  information  from  Figure  (4-1),  but  gives  a  physical  feel 
for  the  parameter  space.  The  stability  transition  curves  obtained  in  Chapter  5  should  be 
plotted  along  this  surface  in  3 -space. 

FLOQUET  THEORY  AND  STABILITY  TRANSITIONS  CURVES 

Equation  (4-15),  the  variational  equation,  is  a  Hill's  equation,  meaning  it  is  a  linear 
differential  equation  with  periodic  coefficients.  This  has  two  important  consequences  in 
relation  to  examining  the  stability  of  the  solutions  of  Equations  (4-1)  along  the  ZMD 
surface. 


Since  Equation  (4-15)  is  a  Hill's  equation  whose  periodic  coefficients  have  period  n,  all 
stable,  periodic  solutions  will  have  period  n  or  2k.  Recalling  that  u(t)  is  periodic,  with 
period  2%,  u2(t)  must  have  period  n,  then  the  claim  follows  directly  from  Floquet's 
Theorem  [16]. 
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The  second  result  of  the  variational  equation  being  a  Hill's  equation,  is  that  the  transition 
from  stable  to  unstable  behavior  in  the  parameter  space  will  always  be  through  a  stable, 
periodic  solution.  In  combination  with  the  first  result,  that  stable,  periodic  solution  must 
be  of  period  tu  or  2n.  To  prove  this  claim,  consider  Floquet's  Theorem[14]: 

The  regular  system  x  =  P(t)x,  where  P  is  an  n  x  n  matrix  function  with 
minimal  period  T,  has  at  least  one  non-trivial  solution  x  =  x(t)  such  that 
x(t  +  T)  =  ux(t)  where  u  is  a  constant. 

The  variational  equation  can  be  written  in  the  form  required  by  Floquet's  Theorem  by 
selecting: 


(4-20) 

and  identifying: 

0  1 

(4-21) 


0  1 

(l  +  2eA)    8(1-28^  -u2) 


which  has  period  n  due  to  the  presence  of  u2(t).  Therefore  Equation  (4-15)  meets  the 
conditions  of  Floquet's  Theorem.  Now  consider  the  possible  values  for  u..  If  u>l  the 
solution  x(t)  becomes  unbounded  over  time  since  the  solution  is  multiplied  by  [i  for  each 
interval  of  time  k,  similarly  if  n<- 1  the  solution  will  also  be  unbounded  over  time.  In  both 
of  these  cases  the  solutions  to  the  variational  equation  are  unstable.  If  -1<u<1,  the 
solution  is  stable  and  decays  to  the  origin  over  time.  That  leaves  two  possible  values  for 
\x.  \i=\  or  u=-l. 

If  u=l,  the  solution  is  periodic  with  period  k,  the  same  period  as  P(t).  If  u=-l,  the 
solution  is  also  periodic  but  has  period  2tc.  Identifying  where  the  characteristic  number,  \i, 
takes  on  the  values  of  1  and  -1  marks  the  transition  between  stable  and  unstable  solutions 
of  the  variational  equation. 

The  fact  that  the  transition  between  stable  and  unstable  behavior  is  always  through  a 
stable,  periodic  solution  with  period  7i  or  2n  focuses  the  study  of  this  thesis.  Locating  the 
periodic  solutions  of  Equation  (4-15)  on  the  ZMD  surface  identifies  the  stability  transition 
curves,  where  the  stability  transition  curves  are  the  curves  which  divide  stable  and  unstable 
behavior  of  Equations  (4-1)  in  the  parameter  space  along  the  ZMD  surface.     Hence,  the 
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next  task  is  to  identify  where  periodic  solutions  (of  period  k  or  2%)  of  Equation  (4-15) 
exist  on  the  ZMD  surface. 


CHAPTER  5:  PERIODIC  SOLUTIONS  OF  THE  VARIATIONAL  EQUATION 


As  discussed  at  the  end  of  Chapter  4,  identifying  the  periodic  solutions  of  the  variational 
equation  along  the  Zero  Mean  Damping  (ZMD)  surface  is  equivalent  to  finding  the 
stability  transition  curves  for  the  coupled  oscillators  modeled  by  Equations  (4-1).  The 
next  task  is  to  find  those  periodic  solutions 

SOLUTION  METHOD 

The  periodic  solutions  of  the  variational  equation  will  be  identified  using  Lindstedt's 
method  [14],  which  was  also  used  in  Chapter  1  to  solve  the  van  der  Pol  equation.  The 
variational  equation  along  the  zero  mean  damping  surface  is  restated  here  for  easy 
referencing: 

v-e(l-2B-u2)v  +  (l  +  2£A)v  =  0  (5-1). 

Some  characteristics  of  the  Equation  (5-1)  are  worth  mentioning.  The  most  significant 
feature  is  that  unlike  the  van  der  Pol  equation  itself,  the  variational  equation  is  a  linear 
ordinary  differential  equation.  Since  u(x)  can  be  determined  to  a  suitably  large  order  in  e 
as  outlined  in  Chapter  1,  and  BZmd  is  the  expansion  for  the  ZMD  surface,  all  the 
coefficients  in  the  variational  equation  are  known  except  8  A.  As  a  consequence  of  being  a 
linear  equation,  solutions  to  Equation  (5-1)  will  also  obey  the  principle  of  superposition. 

Even  though  Equation  (5-1)  is  a  linear  differential  equation,  Lindstedt's  method  is  used  to 
find  v(x)  because  of  the  presence  of  u(x)  which  is  a  power  series.  In  keeping  with 
Lindstedt's  method,  the  first  step  is  to  rescale  time  such  that  x  =  co(e)  t.  Substituting  this 
rescaling  into  Equation  (5-1)  and  changing  the  differentiation  to  the  new  time  scale: 

v"  -  e(l  -  2B7MD  -  u2  )ov'  +  co 2  (l  +  2sA) v  =  0  (5-2) 

where  primes  represent  differentiation  with  respect  to  the  rescaled  time  t.  Notice  that  here 
the  function  co(e)  was  already  determined  in  solving  Equation  (1-1)  and  is  given  by 
Equation  (1-19).  The  frequency  expansion  must  be  the  same  as  for  the  original  van  der 
Pol  solution  because  the  solutions  to  Equation  (5-1)  must  have  the  same  period  or  twice 
the  period  of  its  coefficients;  this  was  discussed  at  the  end  of  Chapter  4.  Plugging  in 
power  series  expansions  for  ©,  v,  BZmd,  u  and  eA  according  to: 
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n 

eA  =  A0+£Aisi 


n.2 


BZMD   ~B0  +2-BiE' 
i=2 

u(t)  =  u0(t)  +  2ui(t)Ei  (5-3) 

v(t)=v0(t)+i;vi(t)8i 

i=l 
n,2 

CO  =  l  +  ^COjE1 


and  collecting  powers  of  e  gives  a  series  of  linear,  ordinary  differential  equations,  the  first 
three  of  which  are  presented  here: 

v;'  +  (l  +  2A0)v0  =  0 

v{'+(l  +  2A0)Vl  =  -2A, v0  +  v;  - 2B0v;  -  u2,  vj  (5-4). 

v2'  +  (l  +  2A0)v2  -  -2A2v0  -2AlV]  -2u0u,v;  +  v{  -2B0v{  -uj  w[  -2co2u;' 

Notice  that  the  left  hand  side  of  each  of  Equations  (5-4)  is  of  the  same  form,  namely: 
LHS  =  v1"+(l  +  2A0)vi  (5-5). 

Insisting  that  the  solution  of  Equation  (5-1)  be  periodic  with  period  %  or  2rc,  as  discussed 
above,  the  coefficient  for  Vi  in  Equation  (5-5)  must  be  an  integer  squared.  This  follows 
directly  from  the  Floquet  theory  discussed  at  the  end  of  Chapter  4,  the  solution  to 
Equation  (5-5)  is  Co^t^/l  -  2A0 )  which  can  produce  v(t)  with  period  n  or  2n  only  if 

A/l-2A0  is  an  integer,  therefore: 

k2  -1 
l  +  2A0  =  k2     or    A0=— — -    where    k  e{l,2,3,...}  (5-6). 

At 

This  gives  rise  to  an  infinite  set  of  starting  values  for  Ao,  and  possible  stability  transition 
curves  out  of  each  starting  value.  The  first  few  Ao's  are  given  here: 

A0  =  U,0,f,4,f,12,...}  (5-7). 
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DEALING  WITH  Ao=-l/2 

When  Ao  =  -1/2,  the  variational  equation  changes  form  significantly,  the  zeroeth  order 
equation  of  Equations  (5-4)  becomes: 

v?  =  0  (5-8) 

which  has  solution: 

v0=C0+C01t  (5-9) 

However,  since  the  goal  is  to  find  periodic,  stable  solutions  to  the  variational  equation,  Coi 
=  0  so  v0  =  Co  is  the  solution.  Plugging  in  this  v0  and  examining  the  first  order  equation 
yields: 

v1»=-2AICfl  (5-10). 

In  order  for  Vi  to  be  bounded,  either  Co  =  0  or  Ai  =  0.  The  former  yields  a  trivial  order 
zero  solution,  but  the  latter  is  a  valid  method  of  solving  Equation  (5-8)  with  a  bounded 
solution.  By  similar  arguments  used  about  Equation  (5-9)  it  can  be  seen  that  vi  =  Ci. 
Following  this  pattern,  the  same  equation  arises  at  each  order  of  e,  so  that  the  stability 
transition  curve  from  Ao=-l/2  is  in  fact  just  eA=  -1/2. 

SOLVING  FOR  THE  STABILITY  TRANSITION  CURVES 

For  the  remainder  of  the  possible  values  for  Ao  listed  in  Equation  (5-7),  the  solution 
procedure  is  a  straight-forward  application  of  Lindstedt's  method.  Whereas  the  solution 
of  the  van  der  Pol  equation  in  Chapter  1  resulted  in  fixing  the  frequency  of  the  response  as 
a  function  of  e,  solving  the  variational  equation  for  v(t)  will  give  the  coefficients  for  the 
eA  expansion. 

Unfortunately,  unlike  the  arguments  made  in  Chapter  1  for  the  van  der  Pol  equation,  there 
are  no  obvious  symmetry  arguments  to  use  in  the  solution  of  the  variational  equation. 
While  the  solution  of  the  variational  equation  shows  interesting  patterns  for  some  values 
of  Ao,  the  overall  pattern  is  not  easily  discernible.  This  results  in  a  very  costly 
computational  procedure  for  determining  the  stability  transition  curves. 

The  stability  transition  curves  calculated  below  and  Tabulated  in  Appendix  E  were 
generated  using  the  Mathematica™  code  found  in  Appendix  B.  One  feature  of  the  code  is 
worth  mentioning  here.  In  an  effort  to  make  solving  the  second  order  differential 
equations  more  efficient,  a  specific  solving  routine  called  MySolve[     ]  was  written.   This 
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routine  takes  advantage  of  the  fact  that  at  every  order  of  8,  the  differential  equation 
encountered  always  takes  the  form: 

v,"+k2Vi  =  f(sin(mt),Cos(mt))  where  m  =  1,  2,  3,...,2i+l  (5-11) 

and  allows  the  code  to  produce  solutions  without  calling  a  general  differential  equation 
solver.  The  specialized  solver  significantly  increases  the  computational  efficiency. 

CONDITIONS  FOR  SUPPRESSING  SECULARITY 

Following  Lindstedt's  method  and  substituting  the  power  series  of  Equation  (5-3)  into 
Equation  (5-2)  and  collecting  terms  of  the  same  order  in  s,  the  first  few  equations  can  be 
obtained: 

k2v[0,  t]+v(0'2)[0,  t]=  0 

k2  v[l,  t]  +v(0'2)[l,  t]  =  -2  a[l]  v[0,  t]  +2  v^fO,  t]  -4Cos[t]2  v(<u)[0,  t] 

k2  v[2,  t]  +  v(0'2)[2,  t]  =  -2  a[2]  v[0,  t]  -  2  a[l]  v[l,  t]  -  3  Cos[t]  Sin[t]  v((U)[0,  t]  + 

Cos[t]  Sin[3  t]  v(01)[0,  t]  +2  v^El,  t]  -4  Cos[t]2  v(<u)[l,  t]  +  -  v(0>2)[0,  t] 

(5-12) 

where  k2  =  1  +  2  Ao  as  defined  in  Equation  (5-6).  The  notation  used  in  Equations  (5-12) 
is  consistent  with  the  Mathematica™  code  found  in  Appendix  B.  In  this  notation,  v[i,t]  is 
the  same  as  vj(t),  a[i]  is  A;,  and  the  numbers  inside  parenthesis  represent  differentiation 
with  respect  to  the  variable  t. 

Noting  that  the  left  hand  side  of  each  of  Equations  (5-12)  are  identical,  it  can  be  seen  that 
the  homogeneous  solution  for  each  order  of  e  will  be  similar: 

vh[i,  t]  =  c[i]  Cos[kt]  +  s[i]  Sin[kt]  (5-13) 

where  c[i]  and  s[i]  will  be  used  consistently  to  represent  the  coefficients  of  the 
homogeneous  solutions  of  the  ith  order  equation.  The  particular  solutions  will  not  be  as 
consistent,  but  will  be  of  the  form: 

2i+l 

vp[i,  t]  =  ^  ( c[i,  m]  Cos[m  k  t]  +  s[i,  m]  Sin[m  k  t] } 

"*>  (5-14). 

The  c[i,m]  and  s[i,m]  are  tabulated  in  Appendix  D. 
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In  the  course  of  solving  Equations  (5-1 1),  two  unique  solutions  for  suppressing  secularity 
will  arise  for  each  Ao  [20].  The  outcome  of  these  two  possible  solutions  for  v[t]  and  8  A  is 
that  for  each  starting  value  of  Ao,  other  than  -1/2,  there  will  be  two  curves  along  the  ZMD 
surface  which  have  valid,  periodic  solutions,  with  period  k  or  2k.  These  curves  intersect 
at  eA  =  Ao  and  resemble  a  horn  with  its  point  downward  and  will  sometimes  be  referred 
to  as  the  "stability  horn"  arising  from  each  Ao.  Extending  this  result  into  the  three 
dimensional  parameter  space  of  Figure  (4-2),  the  two  dimensional  surface  where  periodic 
solutions  to  Equation  (5-2)  exist  will  resemble  cones  will  their  points  coming  out  of  each 
Ao.  The  regions  within  the  cones  are  unstable,  that  is  the  solutions  to  the  variational 
equation  will  grow  over  time,  making  the  origin  of  the  variational  space  an  unstable 
equilibrium  point.  The  regions  outside  of  the  cones  and  above  the  ZMD  surface  will  be 
asymptotically  stable;  all  solutions  will  decay  towards  the  origin  of  the  variational  space 
over  time.  The  actual  surfaces  of  the  cones  are  configurations  which  yield  orbitally  stable 
solutions. 

One  of  the  main  difficulties  in  solving  Equations  (5-12)  will  be  to  identify  this  split  in  the 
stability  conditions.  The  conditions  are  different  for  each  Ao  and  will  be  considered  one  at 
a  time. 

Case  1 :  Ao  =  0 

For  the  case  of  Ao  =  0,  Equation  (5-6)  yields  k=l.  Solving  Equations  (5-12)  and 
suppressing  the  secular  terms  appearing  on  the  right  hand  side  of  the  first  order  term 
gives: 

-c[0]-2a[l]s[0]=0 

-2a[l]c[0]-s[0]=0  (5_15) 

which  has  two  possible  solutions  for  a[l]  and  s[0]: 

{{s[0]  ->  -c[0],  a[l]  -»  -},  {s[0]  -*  c[0],  a[l]  ->  --}) 

U  2J    l  2>>  (5-16). 

The  first  solution  will  yield  a  curve  to  the  right  of  the  second  when  viewed  on  the  ZMD 
surface.  From  this  point  on,  the  solutions  for  v[t]  and  eA  are  unique  and  the  stability 
transition  curves  are  given  by: 

eAow=-y+  h~ ~  ^T"T7TT-- " +  „.. +  0(0 


e2   e3 

+  — 

8   32 

7  e4   3  e5 
384   1024 

247  c6    11657e7    1224811c8 
442368   21233664   2548039680 

£   e3 

8   32 

7c4   3e5 
384   1024 

247  c6    11657e7    1224811c8 
442368   21233664   2548039680 

(5-17) 

6 

:     —  +"  +  -  +  ~  ~+     0(0 
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The  power  series  given  by  Equation  (5-17)  were  calculated  out  to  0(e33)  and  converted 
into  Pade  approximants  which  are  plotted  in  Figure  (5-1)  with  the  power  series. 


-l 


Pade  Approximants 


Power  Series 


Figure  (5-1)  Power  series  and  Pade  approximants 
[16/16]  to  transition  curves  out  of  Ao=0. 


Case  2:  Ao  -  3/2 

For  Ao  =  3/2,  Equation  (5-6)  gives  k=2.    Suppressing  the  secular  terms  in  the  first  order 
equation  of  Equation  (5-12)  yields: 


-  2  a[l]  s[0]  =0 
-2  a[l]  c[0]  =  0 


(5-18) 


which  requires  a[l]=0  for  a  non-trivial  solution.     Suppressing  secularity  in  the  second 
order  equation: 


2  s[0] 
3 
c[0] 


-  2  a[2]  s[0]  =  0 
-2  a[2]  c[0]  =0 


which  has  three  possible  solutions: 

|{a[2]  -►  --,  s[0]  ->  u},  |a[2]  -►  -,  c[0]  -►  o},  {s[0]  ->  0,  c[0]  ->  0}} 


(5-19) 


(5-20). 
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The  last  solution  is  the  trivial  solution,  but  the  first  two  again  provide  the  split  for  the 
lower  and  upper  stability  curves.  From  this  point,  the  solutions  are  unique  and  the  stability 
transition  curves  are  given  by: 


3       e2       109  e4 

cAlow  = + 

2       6        3456 

53  c4 


3419  c6         297305  e8 


10, 


1990656       573308928 


+  0(e,u) 


^up  =  T+_ 


6983^    740213  e8      10 
+  -       -+0(e10) 


3456   1990656   2866544640 


(5-21). 


The  Pade  approximants  for  the  power  series  of  Equations  (5-21)  calculated  to  0(e33)  are 
plotted  in  Figure  (5-2)  with  the  power  series. 
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Figure  (5-2)  Power  series  and  Pade  approximants 
[12/12]  to  transition  curves  out  of  Ao=3/2. 


Case  3:  Ao  =  4 


For  Ao  =  4,  Equation  (5-6)  gives  k=3.    Suppressing  the  secular  terms  in  the  first  order 
equation  of  Equation  (5-12)  yields: 


-2  a[l]  s[0]  =0 
-2a[l]  c[0]  =  0 


(5-22) 


which  requires  a[l]=0  for  a  non-trivial  solution.     Suppressing  secularity  in  the  second 
order  equation: 
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9s[0] 

— —  -  2  a[2]  slO]  =  0 
16 

9c[0] 


-2a[2]  c[0]=0 
16  (5-23) 

which  requires  a[2]=-9/32  for  a  non-trivial  solution.    Suppressing  secularity  in  the  third 
order  equation: 

5  c[0] 

L-L-2  a[3]  s[0]  =  0 

32 

5  s[0] 
-2a[3]  c[0]+  — ±-Uo 

32  (5-24) 

which  has  two  possible  solutions  for  a[3]  and  s[0]: 

{{s[0]  -*  -c[0],  a[3]  ->  -— },  {s[0]  ->  c[0],  a[3]  -*  -^}} 


(5-25). 


These  two  solutions  provide  the  split  for  the  lower  and  upper  stability  curves.   From  this 
point,  the  solutions  are  unique  and  the  stability  transition  curves  are  given  by: 

9  e2       5  c3       2877  e4 

^Aow  =  4 +  + 

32         64        40960 

2537  c5        54707  c6         11663959  c7         17260922723  c8  g 

+ +  O  (<T) 

737280       15728640       14155776000       12683575296000 

9  c2       5  c3       2877  c4       2537  c5        54707  c6         11663959  c7 
cAlin  =  4 + + 


*  32         64        40960        737280       15728640       14155776000 

17260922723  c8  g 

12683575296000  (5-26). 

The  Pade  approximants  for  the  power  series  of  Equations  (5-26)  calculated  to  0(s33)  are 
plotted  in  Figure  (5-3)  with  the  power  series. 
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Figure  (5-3)  Power  series  and  Pade  approximants 
[16/16]  to  transition  curves  out  of  Ao=4. 


CHAPTER  6:  SUMMARY  AND  FUTURE  WORK 


In  Chapter  5,  the  first  three  pairs  of  stability  transition  curves  along  the  Zero  Mean 
Damping  (ZMD)  surface  were  determined.  Combining  Figures  (5-1,2  and  3)  yields  the 
stability  regions  for  the  smaller  values  of  e A.  In  Figure  (6- 1 ),  U  marks  regions  that  are 
unstable  and  S  marks  regions  where  orbitally  stable  solutions  to  the  variational  equation, 
Equation  (4-15),  exist.  The  stability  transition  curves  are  approximated  by  Pade 
approximants. 


Figure  (6-1)  Stability  regions  for  coupled 
VDP  oscillators  on  the  ZMD  surface. 

The  three  dimensional  version  of  the  same  graph  is  given  in  Figure  (6-2)  where  the  three 
dimensions  are  the  three  parameters  of  Equation  (4-1),  {e,  sA,  eB}  .  This  is  a  graph  of  the 
stability  transition  curves  found  in  Chapter  5  plotted  on  the  ZMD  surface  of  Figure  (4-2). 

SUMMARY 


Using  the  techniques  outlined  in  Chapter  5  and  employing  a  symbolic  mathematics 
package  (Mathematica™)  the  first  three  pairs  of  stability  transition  curves  for  a  pair  of 
weakly  coupled  van  der  Pol  oscillators  were  determined.  All  calculations  were  carried  out 
exactly,  that  is  exact  fractions  were  used  in  all  computations  with  no  numerical  roundoff. 
Using  this  analytical  approach,  the  six  stability  transition  curves  seen  in  Figures  (6-2)  and 
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(6-2)  were  calculated  to  0(e33).  The  power  series  coefficients  for  the  transition  curves  are 
tabulated  in  Appendix  E. 

FUTURE  WORK 

There  are  three  main  areas  for  future  work  on  this  project. 

The  first  is  to  extend  the  range  of  eA  covered  in  this  thesis,  namely  find  the  stability 
transition  curves  coming  from  the  rest  of  the  Ao  values  given  in  Equation  (5-7).  Work  on 
this  area,  using  the  code  written  for  this  thesis,  is  already  in  progress.  The  main  difficulty 
in  finding  the  stability  transition  curves  coming  from  the  next  few  values  of  Ao  is  showing 
that  the  curves  are  unique.  Present  conjecture  is  that  there  are  more  conditions  needed  to 
fully  specify  the  upper  curves  than  used  in  Chapter  5.  To  find  what  those  conditions  might 
be,  research  is  underway  to  use  Hill's  determinants  to  find  the  periodic  solutions  rather 
than  Lindstedt's  method. 

Another  area  for  future  work  was  hinted  at  in  Chapter  4  when  the  ZMD  surface  was 
introduced.  The  full  challenge  for  this  type  of  stability  analysis  would  be  to  find  the 
stability  transition  surfaces  in  the  three  dimensional  parameter  space.  Currently,  numerical 
calculations  of  this  project  are  under  way  by  Storti  and  Reinhall  (REFERENCE>»).  The 
logical  extension  of  this  project  would  be  to  find  the  transition  surfaces  analytically  to 
avoid  numerical  roundoff  errors. 

The  final,  and  ultimate  continuation  of  this  work  is  to  physically  realize  the  coupled 
relaxation  oscillators  modeled  by  Equation  (4-1).  Once  these  oscillators  are  accurately 
predicted  by  the  mathematical  models,  they  could  be  used  to  generate  control  signals  for 
ambulatory  robots,  mechanical  fish  or  similar  mechanisms.  The  goal  of  this  thesis  was  to 
further  understanding  toward  reaching  this  goal. 
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Figure  (6-2)  Three  dimensional  plot  of  the 
stability  transition  curves  along  the  ZMD  surface. 
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APPENDIX  A:  MATHEMATICA™  CODE  FOR  VDP  EQUATION 


This  first  block  of  code  produces  the  first  three  terms  of  the  of  the  VDP  equation. 
Comments  are  included  to  help  clarify  the  procedure.  A  second  set  of  code  is  included 
which  will  start  with  a  certain  number  of  terms  already  in  the  solution  and  calculate  two 
more  terms.  That  code  can  be  used  to  extend  the  solution  as  far  as  the  user  wants  or  until 
the  computer  runs  out  of  memory. 

STARTER  PROGRAM 


$Line  =  <»; 

SetDirectory  [  "c :  \\v*%A\conplex\\cdata"  ] 

(*  pmult  [  f ,  g,  x,  n]   multiplies  power  series  f  by  power  series  g 
where  x  is  the  variable,   and  keeps  only  through  order  n  *) 

n  tvi 
jault[f_,  g_,  x_Synfaol,  n_Integer]   :  =  V  V  3^*3  Coefficient  [  f ,  x,  i]  Coefficient  [g,  x,  j] 

i»0  j=0 

(*  getX[n]    creates  a  power  series  in  e  of  order  n  with  argument  x[i,t]    *) 

n 
getX[n_]  :=  ExpandlJ^xli,  t]  ] 

i»0 

(*  xeven[]    and  xodd[]    assign  x[i,t]   to  be  a  complex  cosine  or  sign  *) 

2i+l 
xseventiJ   :=   Z  A[if  j]  (Z^Z'3) 

A  >2 

2i+l 

xodd[iJ  :=  Z  IA[i'  jl  (z3-2"11) 

A>2 

(*  getW[n]  creates  the  frequency  power  series  *) 
n 

getW[n_j  :=  1+  V  e1  w[i] 

U2 
A  i=2 
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(*  noslopef]  enforces  the  boundary  condition  in  Equation  (1-1)  *) 
noslope[i_]  :=  Module  [{temp} ,  tenp=  Coefficient  [X[t]  ,  e,  i]  ; 
Solve [  (myD[tenp]  /.  Z-»  1)  ==  0,  A[i,  1]]] 

(*  myD[]  and  myD2[]  are  specialized  derivatives  *) 
myD[fj  :=f/.{Z-->IaZa,  Z-»IZ} 
myD2[fj  :=  mvD[myD[f]  ] 


ClearAll[x,  w,  X,  W,  temp,  lhs,  rhs,  A,  B,  i]  ;  (*  start  with  clean  slate  *) 

n=  3;  (*  order  of  solution  *) 

X[t]  =  getX[n]  ;  (*  build  power  series  in  e  *) 
W=  getW[n]  ; 

(*  build  the  right  hand  sides  of  Eq  1-5  *) 

rhs  =  pxnult[e,  pmilt[W,  pnult[dtX[t]  ,  1-  E«ult[X[t]  ,  X[t]  ,  e,  n-  1]  ,  e,  n-  1]  , 

e,  n-1],  e,  n]  ; 
lhs=  Fmult[pmult[W,  W,  e,  n]  -1,  d{t,2>X[t]  ,  e,  n]  ; 

(*  assign  x[i,t]    form  to  match  Eq  1-15  *) 
Table[x[i,  t]  =  xeven[i]  ,  {i,  0,  n,  2}]  ; 
Table[x[i,  t]  =  xodd[i] ,  {i,  1,  n,  2}]  ; 

(*  back  to  building  right  hand  sides  *) 

ndd=  rhs-  lhs  /.  Table[x(0,1>  [i,  t]  -♦  ra^D[x[i,  t]]  ,  {i,  0,  n}]  ; 
mid=  mid/.  Table  [  x(0 ' 2)  fi,  t]  -»myD2[x[i,  t]  ]  ,  {i,  0,  n}]  ; 
X[t]  =  X[t]  ; 

(*  pull  out  each  right  hand  side  for  each  order  of  e  *) 
tenp=  Table  [Coefficient  [mid,  e,  i]  ,  {i,  1,  n}]  ; 

(*  free  up  some  memory  *) 
CleaxAll [ lhs ,  rhs,  mid,  x]  ; 
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( *  Solving  the  1  st  order  equation  * ) 

i=  l; 

right  =  Collect[ Exparxif  ^^i1!  ]  ,  Table[Z2j+1,  {j,  0,  i}]]  ; 

sec=  E3qjand[  Coefficient  [right,  Z]  ]  ;     (*   Identifies  the  secular  terms  *) 

solRule  =  Solve  [  sec  =  =  0,  A[0,  1]]  ;  (*  Suppresses  secularity  *) 

A[0,  1]  =  A[0,  1]  /.  solRulej3j  ; 

right  =  right  /  .  solRule|[3]  ; 

Coefficient  [right,  Z23+1] 
Table[A[i,  2  j  +  1]  = 2 '  tD'     '  XJ1;    (*   E^01063  E<3  1_15  *> 

A[i,  1]  =  A[i,  1]  /.  noslope[i]  |[1]| ;  (*  Enforces  boundary  condition  *) 

tenplij  =  0; 


( *   Solving  the  2  nd  order  equation  * ) 

i=2; 

right  =  Collect[Expand[tenp|[iJ]  ,  Table[Z23+1,  {j,  0,  i}]]  ; 

solRule=  Solve  [Expand[  Coefficient  [right,  Z]  ]  ==  0,  w[i]]  ;  (*  Id'  s  the  secular  tenns  *; 

w[i]  =  w[i]  /  .  solRuleffl]  ; 

right  =  right  /  .  solRulelll  ; 

Coefficient  fridht,  Z23+1] 

TablefA[i,  2  j  +  1]  = — =-— ,  {j,  1,  i}];  (*  Enforces  Eq  1-15  *) 

1  -  (2  j+  l)2 

tenpHi]  =  0; 


( *  Solving  the  3  rd  order  equation  * ) 
i=3; 

right  =  CollectfE^andt^^1^-]  ,  Table[Z23+1,  { j,  0,  i}]]  ; 

solEule=  Solve  [Expand[  Coefficient  [right,  Z]  ]  ==  0,  A[i-  1,  1]]  ;  (*  Id'  s  the  secular  terms  *) 

A[i-  1,  1]  =A[i-  1,  1]  /.  solRuleElJ; 

right  =  right  /  .  solKule[l]  ; 

,   r   .         Coef ficient[ right,  Z23*1]        .  ,     „  _      „  -  -- 

Table  A[i,  2  j  +  1]  = ^— ,  {j,  1,  1}  ;  (*  Enforces  Eq  1-15  *) 

1-  (2  j+  l)2 

A[i,  1]  =  A[i,  1]  /.  noslope[i]  [1]  ;  (*  Enforces  boundary  condition  *) 

tenpHi]  =  0; 
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( *  Saving  the  data ...  * ) 
DeleteFile  [  "af ile .  txt"  ] 
Save  [  "af  ile .  txt" ,  A] 
DeleteFile  [ "  wf  ile .  txt"  ] 
Save ["wfile.txt",  w] 
DeleteFile  [  "nf ile .  txt"  ] 
Save[  "nf ile .  txt" ,  n] 

(*  This  section  puts  the  coefficients  back  in  a 

form  to  use  with  sines  and  cosines . . . * ) 
Table[Table[B[i,  j]  =  2A[i,  j]  ,  {j,  1,  2x+l,  2}],  {i,  0,  n,  2>]  ; 
Table[Table[B[i,  j]  =  -2A[i,  j]  ,  {j,  1,  2i+l,  2}],  {i,  1,  n,  2}]  ; 
Save[ "bfile.txt",  B] 

(*  Rebuilding  X[t]  in  terms  of  sines  and  cosines  *) 

2i+l 

Table[x[i,  t]  =  £  Bl±r  J]  Cos[j  t]  ,  {i,  0,  n,  2}]  ; 

3»1 
A  j=2 

21+1 

Table[x[i,  t]  =  £  B[i,  j]  Sin[jt] ,  {i,  1,  n,  2}]  ; 

j=i 

A  3=2 

X[t]  =  Sun[x[i,  t]  ex,  {i,  0,  n}]  ; 

(*  Check  out  the  work...   *) 

W 

X[t] 


CONTINUATION  PROGRAM 


The  continuation  program  uses  the  same  functions  defined  above,  and  the  code  is  identical 
down  through  defining  noslope|  ]  then  change  the  code  as  follows: 
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ClearAll[x,  w,  X,  W,  tenp,  lhs,  rhs,  A,  i,  n]  ; 
n«  "nfile.txt"; 
old=  n; 
n  =  old+  2  ; 

Print  ["Starting,  old  =   ",  old,  "     new  =   ",  n] 
tl  =  TiiteUsed[  ]  ; 
X[t]  =  getX[n]  ; 
W=  getW[n]  ; 

rhs  =  l-rnult[X[t] ,  X[t]  ,  e,  n-1]  ; 
tenp=  pnult[W,  dtX[t]  ,  e,  n-  1]  ; 
rhs  =  pmlt[rhs,  tenp,  e,  n  -  1]  ; 
Clear  [tenp]  ; 
rhs  =  rhs  /  .  Tablefe1  -»  a1 ,  {i ,  old,  n  -  1}  ]  ; 

rhs  =  rhs/  .  e-»  0; 
rhs=  rhs  /  .  a-»  e; 
rhs  =  pnult[e,  rhs,  e,  n] ; 
lhs=  pnult[W,  W,  e,  n]  ; 
lhs=  lhs-1; 

lhs  =  pnult[lhs,  d{t,2}  x[t]  t  e/  nl  ! 

lhs=  lhs/.  Tablele^^^a1,  {i,  old+  1,  n}]  ; 

lhs=  lhs/.  e-»  0; 
lhs  =  lhs  / .  a-»  e; 
Table[x[i,  t]  =  xeven[i]  ,  {i,  0,  n,  2}]  ; 

Table[x[i,  t]  =  xodd[i] ,  {i,  1,  n,  2}]  ; 

mid=  rhs- lhs  /.  Table[x(0,1)  [i,  t]  -♦  myD[x[i,  t]  ]  ,  {i,  0,  n}]  ; 

mid=  mid/.  Table[x(0,2)  [i,  t]  -»myD2[x[i,  t]]  ,  {i,  0,  n}]  ; 

Print[ "v?)dating  X[t]  . . ."] 

X[t]  =  X[t]  ; 

Print [  "building  tenp "] 

tenp=  {mid/  .  {e01^1  -*  1,  e-  0}  ,  mid/  .  {e01*2-  1,  e-  0}}  ; 
ClearAll[lhs,  rhs,  mid,  x]  ; 
A«  "afile.txt"; 
w  <<  "wfile.txt"; 
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Tahle[A[i,  2  j  +  1]  = n      ,oi.i2 '  {j'  1'i>l; 


i  =  old+  1  ; 

rhs  =  Collect [Expand[tenpl  11]  ,  TablefZ23*1,  { j,  0,  i}]]  ; 

solRule  =  Solve [ Expand [ Coefficient [rhs,  Z]  ]  ==  0,  w[i]]  ; 

w[i]  =  w[i]  /  .  solRule|[l]|  ; 

rhs  a  rhs  /  .  solRulejl]]  ; 

Coefficient  [rhs,  Z2j+1] 
l-(2j+l): 
tenp[l]|  =  0; 
i  =  old+  2  ; 

rhs =  Collect [Expand[  tMpl  J  ]  ,  Table[Z2j+1,  {j,  0,  i}]]  ; 

solRule=  Solve  [Expand[  Coefficient  [rhs,  Z]  ]  ==  0,  A[i-  1,  1]]  ; 
A[i-  1,  1]  =A[i-  1,  1]  /.  solRulejl); 
rhs  =  rhs  /  .  solRule|[l]  ; 

Coefficient  [rhs,  Z23+1] 
l-(2j+l): 
A[i,  1]  =A[i,  1]  /.  noslope[i]|[ll|; 
tenp|[2l  =  0; 

DeleteFile[{  "afile.txt",  "wfile.txt",  "nfile.txt"}] 
Save  [  "af ile .  txt" ,  A] 
Save  ["wfile.txt",  w] 
Save[ "nfile.txt",  n] 


From  here  the  last  few  routines  in  the  starter  program  can  again  be  used  to  translate  the 
solution  back  into  sines  and  cosines  and  save  the  results  in  that  form. 


Table[A[r,  2j  +  l]  = ^     /2 l-  ,  {j,  1,  i}]  ; 


APPENDIX  B:  STABILITY  TRANSITION  CODE 


The  first  program  is  "Recover  nb"  which  converts  the  coefficients  of  the  complex  VDP 
solution  into  real  coefficients  for  use  in  the  stability  routines.  It  also  calculates  the  ZMD 
surface  coefficients.  It  updates  these  coefficients  in  the  c:\stabil\data  folder  as  well  as  the 
frequency  coefficients.  This  code  needs  to  be  run  only  once  after  generating  as  many 
terms  in  the  VDP  solution  as  necessary. 


RECOVERNB 

Directory  structure  is  assumed  to  be 
c: 

-complex  THIS  HAS  THE  VDP  SOLUTION  STORED  IN  IT 

-  stabil  DATA  FOR  THE  STABILITY  CURVES 

$Line=  a>; 

SetDixectory  [  "c :  \\vcfc>"  ]  ; 

n<<  "ccnplex\\adata\\nfile.txt";         (*  nutter  of  terms  available  in  VDP  solution  *) 

Reads  in  solution  to  complex  solution  and  converts  them  to  coefficients  for  real  solution 

A  <  <  "oarplex\\odata\\af ile .  txt" ; 

Table[Table[B[i,  j]  =  2A[i,  j]  ,  {j,  1,  2i  +  l,  2}]  ,  <i,  0,  n,  2}]  ; 

Tablet  Table  [B[i,  j]  =  -2A[i,  j]  ,  {j,  1,  2i+l,  2}],  {i,  1,  n,  2}]; 

DeleteFile  [  "oaqplex\\adata\\bf ile .  txt"  ]  ; 

DeletePile["stabil\\data\\bfile.txt"]  ; 

Save t "oacplex\\odata\\bf ile. txt",  B]  ; 

Save[ "stabilWdataWbfile.txt",  B]  ; 

The  next  four  lines  calculate  the  coefficients  of  the  ZMD  surface 

2i+l  1 

Table[x[i]  =  ^  -  Bti,  Jl  (Zj  +  Z'3)  ,  {i,  0,  n,  2}]  ; 

>i    2 

A>2 

2i+l 


Table[x[i]  =  2  "-  IBU^  J]  (z3  "  z"3)  •  <*■>  1,  n,  2}]  ; 
3»1      2 

A  j=2 

Table[b[i]  =  JJ  (Expand[-—  x[  j]  x[i-  j]  ]  /  .  Z— »  0)  ,  {i,  0,  n,  2}]  ; 


J-0 

b[0]  =b[0]  +  — ; 
2 


Saves  the  ZMD  coefficients,  updates  the  frequency  coefficients  in  the  stabil  folder  and 
updates  how  many  terms  are  available  in  the  VDP  for  the  stability  routines 
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DeleteFile  [  "stabil\\data\\crit .  txt"  ] 
Save[  "stabil\\data\\crit .  txt" ,  b] 
DeleteFile  [  "stabil\\data\\wf  ile .  txt"] 

CcpyFile [  "ccnplex\\odata\\wf ile .  txt" ,  "stabil\\data\\ wf ile .  txt" ] 
DeleteFile  [  "stabil\\data\\nf  ile .  txt"] 
Save[  "stabil\\data\\nf ile .  txt" ,  n] 

The  next  program  will  start  the  solution  for  the  stability  transition  curve: 


STARTERNB 

Assumes  the  same  directory  structure  as  Recover,  nb 

$Line=  oo 

SetDirectory[  "c :  \\vdp\\stabil\\data"] 

pmult[  ]  is  the  same  as  in  the  VDP  code.   getV[  ]  is  equivalent  to  getX[  ]  from  VDP  and 
getA[  ]  makes  the  power  series  for  the  transition  curve 

n  n-i 

pnult[f_,  g_,  xJSpibol,  n_Integer]   :=  V  Vx1*3  Coefficient[f ,  x,  i]  Coefficientfg,  x,  j] 

i.0>0 


getVforderj  :=    V  env[n,  t] 
n=0 
crac 
getA[order_j  :=    V  ena[n] 

nrO 

update[  ]  updates  the  RHSs  of  Equation  (5-12)  as  each  v[i,t]  is  determined. 

v?riate[tenp_,  i_]  :=  Module[{ttt} ,  ttt=  tenp; 
tttli-l]  =0; 

Table[ttt|[il  =  ttt[il  /.  v[j,  t]  -+  v[j,  t]  ,  {j,  0,  i-  1}]  ; 
Table[ttt[i]|  =  tttlil  /.  v(0'w  [  j,  t]  -etv[j,  t]  ,  {j,  0,  i-1}]  ; 
Tableftttlil  =  tttlil  /.  v(0'2)[j,  t]  ^a{t,2}V[j,  t]  ,  { j,  0,  i - 1}]  ; 
ttt] 

Specialized  ODE  solver  discussed  in  Chapter  5. 
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mysolve[i_,  rightj  :=  Module[{n,  temp} , 

n=  2  (i+  1)  +  1; 

tenp=  c[i]  Cos[kt]  +  s[i]  Sin[kt]  ; 

Z*"1  Coefficient  [right,  Cos[mt]]  Cos[mt] 
m-l  -m2+k2 

A     Coefficient  [  right ,  Cos  [  m  t]  ]  Cos  [  m  t] 

tenp  =  tenp  +     >     — — ; 

i&i  -m2+k2 

^  Coefficient [right,  Sin[mt]]  Sin[mt] 
tenp  =  tenp+   > 


-m2+k2 


A     Coefficient t  rigjit ,  Sin[mt]  ]  Sin [mt] 
tenp  =  tenp+     > 


•Si  -m2+k2 


tenp] 


Build  the  power  series  for  the  frequency,  W,  the  VDP  solution,  U,  and  the  ZMD  surface, 
B 

5 

W=  1+    Ye1  w[i]  ; 


i=2 
A  1.2 

5 


U=2eiu[i,  t]; 

i=0 
5 

B=  £  e^bii]; 


1=0 

ClearAll[n,  v,  V,  A,  a,  tenp,  right,  s,  c] 

n=  5;  <*  how  many  terms  to  calculate  *) 

V[t]  =  getV[n]  ; 

A=  getA[n]  ; 

a[0]  =  0;  (*  which  value  of  a[0]   to  go  after  *) 


k=Vl+2a[0]  (*  k  of  Eq  (5-6)    *) 

rhs  is  the  right  hand  side  of  Equation  (5-12)  prior  to  separating  the  powers  of  e. 

Once  they  are  separated,  each  power  is  stored  in  temp[i]. 

rhs=  pmult[pmult[eW,  atV[t]  ,  e,  n]  ,  1  -  pmult[U,  U,  e,  n-  1]  -  2B,  e,  n]  - 

2pmult[A-aO,  V[t]  ,  e,  n]  -pmult[pmult[W,  W,  e,  n]  -  1,  d(t,2}V[t]  ,  e,  n]  ; 
temps  Table  [Coefficient  [rhs,  e,  i]  ,  {i,  1,  n}]  ; 
Clear [rhs] ; 

Read  in  the  coefficients  for  VDP,  frequency  and  ZMD  surface  (generated  by  recover. nb). 
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u<<  "newu.txt"; 

w<<  "wfile.txt"; 

b«  "crit.txt"; 

Solve  the  zeroeth  order  equation,  note  the  rhs  is  zero  for  this  equation.    Last  line  is  the 
condition  of  unit  displacement  for  the  perturbation.     It  can  be  replaced  with  other 
conditions  or  left  out  entirely. 
i=0; 

v[i,  t]  =  mysolvefi,  0]  ; 

v[i,  t]  =  v[i,  t]  /.  Solve!  (getV[i]  /.  t-»  0)  ==  1,  c[i]]IU; 

Solve  the  first  order  equation.   First  update  the  rhs  to  reflect  v[0,t],  next  find  the  secular 

term  coefficients,  si  and  cl .  Next  solve  to  eliminate  the  secular  terms,  in  this  case  for  a[l] 

and  s[0].  Next  update  the  coefficients  and  the  rhs  with  the  solution  just  found  and  finally 

call  the  specialized  ODE  solver  to  get  v[l,t].  Again  the  last  line  imposes  unit  displacement 

and  could  be  left  out  or  replaced  by  other  conditions. 

Note:  For  this  case,  a[0]=0,  picking  sol[[l]]  gives  the  lower  curve  while  sol[[2]]  would 

give  the  upper  curve. 

i=i; 

tenp=  update  [tenp,  i]  ; 

right  =  Collect [TrigReduoe[tenp[i]|]  ,  {Sin[t]  ,  Cos[t]  }]  ; 

si  =  Coefficient  [right,  Sin[t]  ]  ; 

cl  =  Coefficient  [right,  Cos[t]  ]  ; 

sol=  Solve[{sl==  0,  cl==  0},  <a[i] ,  s[i-l]}] 

a[i]  =  a[i]  /.  sol|[l]|  ; 

v[i-l,  t]  =v[i-l,  t]  /.solilj; 

s[i-l]  =s[i-l]  /.  aol[l]; 

right  ■  right  /  .  sol  [I]  ; 

v[i,  t]  =  mysolve[i,  right]  ; 

v[i,  t]  =  v[i,  t]  /.  Solve[(getV[i]  /.  t->  0)  ==  1,  c[i]]g;lj; 

Solve  second  order  equation,  just  like  above. 

i=2; 

ternp=  update[tenp,  i]  ; 

right  =  Collect [TrigRBduce[tenp|[il]  ,  {Sin[t]  ,  Cos[t]  }]  ; 

si  =  Coefficient  [right,  Sin[t]  ]  ; 

cl  s  Coefficient  [right,  Cos[t]  ]  ; 

sol=  Solve[{sl==  0,  cl==  0} ,  {a[i]  ,  s[i-  1]}] 

a[i]  =  a[i]  /.  solHJ  ; 

s[i-l]  =s[i-l]  /.  sol[ll; 

right  =  right  /  .  sol I IB  ; 
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v[i,  t]  =  mysolvefi,  right]  ; 
v[i,  t]  =  v[i,  t]  /.  Solve [(getV[i]  /  .  t-»  0)  ==  1,  c[i]]  Jlfl  ; 

Same  thing  for  third  order.  . 

i=3; 

tenp  =  update  [tenp,  i]  ; 

right  =  Collect  [  TrigReduce  [tenplifl]  ,  {Sin[t]  ,  Cos[t]}]  ; 

si  =  Coefficient  [right,  Sin[t]  ]  ; 

cl  =  Coefficient  [right,  Cos[t]  ]  ; 

sol=  Solve[  {si  ==  0,  cl==  0> ,  {a[i] ,  s[i-l]}] 

a[i]  =  a[i]  /.  soljlj  ; 

s[i-l]  =s[i-l]  /.  solllj; 

right  =  right  /  .  sol|[l]  ; 

v[i,  t]  =  mysolve [ i ,  right]  ; 

v[i,  t]  =  v[i,  t]  /.  Solve! (getV[i]  /.  t-»  0)  ==  1,  e[i]]IU; 

Same  thing  for  fourth  order. . . 

i=4; 

terrp  =  update!  tenp,  i]  ; 

right  =  Collect [ TrigReduce [ tenpff i| ]  ,  {Sin[t]  ,  Cos[t]}]  ; 

si  =  Coefficient  [right,  Sin[t]  ]  ; 

cl  =  Coefficient  [right,  Cos[t]  ]  ; 

sol=  Solve[{sl==  0,  cl==  0}  ,  {a[i]  ,  s[i-  1]}] 

a[i]  =  a[i]  /.  solIU; 

s[i-l]  =s[i-l]  /.  sollll; 

right  =  right  /  .  sollll  ; 

v[i,  t]  =  mysolve[i,  right]  ; 

v[i,  t]  =  v[i,  t]  /.  Solve[(getV[i]  / .  t->  0)  ==  1,  c[i]]|[l]|; 

Same  thing  for  fifth  order. . . 

i=5; 

tenp=  update  [temp,  i]  ; 

right  =  Collect [ TrigReduce [tenp[ij]  ,  {Sin[t]  ,  Cos[t]  }]  ; 

si  =  Coefficient  [right,  Sin[t]  ]  ; 

cl  =  Coefficient  [right,  Cos[t]  ] ; 

sol=  Solve[ {si  ==  0,  cl==  0},  {a[i] ,  s[i-l]}] 

a[i]  =  a[i]  /.  solilj  ; 

s[i-l]  =s[i-l]  /.  sol[ll; 

right  =  right  /  .  sol|[l]| ; 

v[i,  t]  r  mysolve[i,  right]  ; 

v[i,  t]  =  v[i,  t]  /.  Solve[(getV[i]  /.  t-»  0)  ==  1,  c[i]]  Jlj  ; 
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Save  the  results.  aOl.txt  is  the  lower  curve  out  of  a[0]=0.  nOl.txt  is  how  many  terms  of 
v[i,t]  have  been  solved. 

DeleteFile [ "aOl .  txt"  ] 
Save["a01.txt",  A] 
DeleteFile  [  "vOl .  txt"  ] 
Save["v01.txt",  V) 
DeleteFile [ "nOl . txt" ] 
Save ["nOl. txt",  n] 


CONTINUE  NB 

Continue  will  take  the  results  of  either  starter. nb,  or  itself  and  extend  the  results  by  four 
terms.  The  functions  down  through  and  including  mysolve[  ]  from  starter,  nb  need  to  be 
included  as  well  as  MyExpand[  ]  which  uses  Euler's  identity  to  update  the  right  hand  side 
of  Equation  (5-12)  more  quickly: 

MyExpand[mess_]   :=  Module[{Z,  temp}  , 

tenp=  Expand[mess  / .  {Cos[m_t]  ->  —  (Zm+  Z"m)  ,  Sin[m_t]  -*  —  I  (Zm-  Z'm)  , 

Cosrt]  -,  I  (Z+|)  ,  Sin[t]  -  -|  I  (Z-  I)}]  ; 

Expand[temp  /  .  {Zm--+  Cos[m  t]  +  I  Sin[m  t]  ,  Z2  -»  Cos[2  t]  +  I  Sin[2  t]  , 
Z-+  Cos[t]  +  ISin[t]}]] 
The  next  section  of  code  checks  to  ensure  enough  terms  are  available  in  the  VDP  solution 
and  then  builds  the  right  hand  sides  of  Equation  (5-12)  much  as  above.   It  prints  out  lines 
to  let  the  user  know  what  order  terms  it  is  trying  to  calculate: 

n<<  "nfile.txt"; 

m=  n; 

Print [m,  "  terms  in  VDP  solution  available."] 

n«"n01.txt"; 

min=  n+  1; 

n=  min+  3; 

Print  ["Starting  with  ",  min,  ",  going  to  ",  n,  "."] 

If[n>  m,  Print["Not  enough  terms! !"]  ;Quit[]  ,  Clear [ m]  ] 
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n 

W=  1+   V  e1  w[i]  ; 


1=2 
A  1=2 


U^e^i,  tj; 

i=0 

n 
ZMD=    V  eib[i] 


1=0 

A  1=2 

ClearAll[v,  V,  A,  tenp,  rhs,  s,  c] 

V[t]  =  getV[n]  ; 

A  =  getA[n]  ; 

rhs=  pmult[eW,  atV[t]  ,  e,  n]  ; 

tenp=  1-  pnultfU,  U,  e,  n]  -2ZM3; 

rhs  =  pmult[rhs,  tenp,  e,  n,  min]  ; 

tenp=  2pnult[A- a[0]  ,  V[t]  ,  e,  n,  min]  ; 

tenp2  =  pnult[W,  W,  e,  n]  -  1; 

tenp2  =  pmult[tenp2,  d{t,2)V[t]  ,  e,  n,  min]  ; 

xhs  =  zhs  -  tenp  -  tenp2  ; 

ClearAll[tenp,  temp2]  ; 

tenp=  Table [Coefficient [rhs,  e,  i]  ,  {i,  min,  n}]  ; 

Clear [rhs] ; 

«  "a01.txt"; 

«  "vOl.txt"; 


k=  Vl+  2a[0]  ; 


21+1 

Table[u[i,  t]  =   2  B[i,  j]  Cos[jt]  ,  {i,  0,  n,  2}]  ; 
5-1 

A  3=2 
21+1 

Table[u[i,  t]  =  2  B[i,  3l  Sin[jt],  {i,  1,  n,  2}]; 

3=1 
A  j=2 

B«  "bfile.txt"; 
w<<  "wfile.txt"; 
b«  "crit.txt"; 

The  next  section  calculates  the  four  new  terms,  and  are  just  like  the  equivalent  sections  of 
starter. nb  except  for  using  MyExpand[  ]  instead  of  Expand[  ]: 

i  =  min; 

Print  [  "Working  on  term:    ",  i] 

tenp  =  update  [tenp,  i]  ; 

rhs  =  MyExpand[tenp|[i-  min  +  1]  ]  ; 
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si  =  Coefficient [rhs,  Sin[k  t]  ]  ; 
cl  =  Coefficient  [rhs,  Cos[k  t]  ]  ; 
sol  =  Solve[{sl  ==  0,  cl  ==  0}  ,  {a[i]  ,  s[i  -  1]  }]  ; 
a[i]  =  a[i]  /.  soljlj  ; 
s[i-  1]  =s[i-  1]  /.  solnlj; 
rhs  =  rhs  /  .  sol|[l]|  ; 
v[i,  t]  =  mysolvefi,  rhs]  ; 
v[i,  t]  =  v[i,  t]  /  .  Solve[  (getV[i]  / .  t-»  0)  »  1,  c[i]  ]  ilj  ; 

i=  min+  1; 

Print  ["Vforking  on  term:    ",  i] 

tenp  =  update[teTtp,  i]  ; 

rhs  =  Mty&q3and[  tenp|[i  -  min  +  lj  ]  ; 

si  =  Coefficient  [rhs,  Sin[k  t]  ]  ; 

cl  =  Coefficient  [rhs,  Cos[k  t]  ]  ; 

sol=  Solve[{sl==  0,  cl=  =  0},  {a[i]  ,  s[i-  1]}]  ; 

a[i]  =  a[i]  /.  sol|[l]|  ; 

s[i-  1]  =  s[i-l]  /.  solffU; 

rhs  =  rhs  /  .  solllj  ; 

v[i,  t]  =  mysolve[i,  rhs]  ; 

v[i,  t]  =  v[i,  t]  /  .  Solve[  (getV[i]  / .  t-+  0)  ==  1,  c[i]  ]  [1]  ; 

i=  ndn+  2; 

Print  ["Vforking  on  term:    ",  i] 

tenp=  update  [tenp,  i]  ; 

rhs=  l^xpand[tenp|[i-  min+  1]]  ; 

si =  Coefficient [rhs,  Sin[k  t]  ]  ; 

cl  =  Coefficient  [rhs,  Cos[k  t]  ]  ; 

sol=  Solve[{sl==  0,  cl==  0},  (a[i]  ,  s[i-  1]}]  ; 

a[i]  =a[i]  /.  sol[[l]|  ; 

s[i-  1]  =s[i-l]  /.  sol[[l]|; 

rhs  =  rhs  /  .  sol[lj  ; 

v[i,  t]  =  mysolve[i,  rhs]  ; 

v[i,  t]  =  v[i,  t]  /  .  Solve[  (getV[i]  /  .  t->  0)  ==  1,  c[i]  ]  |[1]|  ; 

i  =  min-t-  3; 

Print [ "Vforking  on  term:    ",  i] 

tenp=  update  [  temp ,  i]  ; 

rhs  =  MyExparid[temp([i-  min-t-  1]]  ; 

si  =  Coefficient  [rhs,  Sin[k  t]  ]  ; 

cl  =  Coefficient  [rhs,  Cos[k  t]  ]  ; 

sol=  Solve[{sl==  0,  cl==  0} ,  (a[i]  ,  s[i-  1]}]  ; 
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a[i]  =  a[i]  /.  solilj  ; 
s[i-l]  =s[i-  1]  /.  solfllfl; 
rhs  =  rhs  /  .  soljlj  ; 
v[i,  t]  =  mysolve[i,  rhs]  ; 
v[i,  t]  =  v[i,  t]  /  .  Solve[  (getV[i]  /  .  t-»  0)  ==  1,  c[i]  ]  [1]|  ; 

Finally  save  the  results  as  above: 

A=  getA[n]  ; 
V[t]  =getV[n]  ; 
DeleteFile  [  "aOl .  txt"  ] 
Sawe["a01.txt",  a] 
DeleteFile  [  "vOl .  txt"  ] 
Save ["vOl. txt",  v] 
DeleteFile  [  "nOl .  txt"  ] 
Save ["nOl. txt",  n] 


APPENDIX  C:  SOLUTION  TO  VDP  EQUATION 


The  analytical  solution  to  Equation  (1-1)  through  order  c'°: 

x[t]  =  2Cos[t]  +e2  (-  Cos[t]   +  _L  Cos[3t]  -  —  Cos [5  t] )  + 
1  8  16  96  ; 

4/73Cos[t]        47Cos[3t]        1085Cos[5t]        2149Cos[7t]        61  Cos  [  9 1] 


e 


12288       1536        27648        110592       20480 
6/6479Cos[t]   48437  Cos[3t]   259945Cos[5 1] 


6635520      35389440       31850496 
4253767  Cos  [  7  t]   480523  Cos  [  9 1]   4937537  Cos  [  11 1]   715247  Cos  [  13 1] 


398131200  73728000        2654208000       3715891200 

377080601  Cos  [t]  36921629  Cos  [  3 1]   4094345  Cos  [  5 1] 

1712282664960      71345111040   +   9172942848 
107409503789  Cos  [  7 1]   188691979247  Cos  [  9 1]   2067847855711  Cos  [11 1] 


45864714240000        59454259200000        936404582400000 
1114925731231  Cos  [  13  tl   784260289  Cos  f  15  tl   392636471  Cos  [17 1]  i 
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1310966415360000       4661213921280      29964946636800   ; 
44898976356847  Cos  [t]  _  4475878408049  Cos  [  3 1] 

1   3020466620989440000     50341110349824000   + 
511229379671  Cos  [  5  tl   21294508407929  Cos  [  7  t]   1838562603094127  Cos  [9 1] 
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2958824445050880      165112971264000000      2621932830720000000 
1299787042060760957  Cos  [  11 1]   11188471201628174701  Cos  [  13 1] 

1321454146682880000000        14800286442848256000000 
27736422934204291  Cos  [  15 1  ]   148525877527816577  Cos  [  17 1] 

78934861028524032000       1522315176978677760000 
1600597132693073  Cos  [  19 1]   29654422883  Cos  [  21 1]  * 

108736798355619840000      32288758824960000   '  + 

,3Sin[t]   1  „.   0^  v   3/  7Sin[t]    21  „.  _     35  „.   r     7   .  _   > 

— Sin[3t]  +eJ —  + Sm[3t] Sin[5t]  + Sin[7tl  + 

1    4     4       I  256    256         576         576       ' 

5  /  12971  Sin[t]   2591Sin[3t] 

'    4423680       294912    + 

52885Sin[5t]   110621  Sin[  7 1]   7457Sin[9t]   5533  Sin[  lit]  \ 


2654208        6635520        1228800       7372800 
.7  /  33114653  Sin[t]   82937Sin[3t]   1939625  Sin[  5 1]   1061235889  Sin[  7 1] 


1   44590694400      212336640       764411904        191102976000 

179467921  Sin[  9 1]   223342933  Sin  [lit]   195050323  Sin[  13 1]   138697  Sin  [15 1] 

+ 


35389440000        92897280000        346816512000       2774532096 
.9  /  9822191352131  Sin[t]   89866136849  Sin[ 3 1]   1702014109  Sin[ 5 1] 

^  251705551749120000  +   342456532992000      12328435187712 
2088126932941  Sin[  7 1]   41070940258801  Sin[  9 1]   1694102483205143  Sin[  11 1] 


2751882854400000       24970788864000000       1048773132288000000 
662514241768639  Sin[  13  tl   507370714969  Sin[  15  tl   15080405867887  Sin[  17  V 


734141192601600000       1740186530611200       302046662098944000 
466445839  Sin[  19 1]  » 

134842259865600  ' 


59 

Frequency  Expansion  through  order  20: 

e2       17  e4         35  e°  678899  e8 

w  =    1 + + 


16   3072   884736   5096079360 
28160413  e10     16729607288111  e12  5722795344767278507  e14 


2293235712000   3698530556313600000   5219366321069752320000000 

1846779765852173498887007  e16      3202506386889338212533493973243  e18 


14731139504587268947968000000000   41577168137747107878744883200000000000 
4038372199485064617691118289043994731  e20 


3872464178615255430139995425341440000000000000 


NUMERICAL  RESULTS 

These  lists  of  numbers  represent  the  coefficients  of  the  cosine  or  sines  of  the  VDP 
solution.  c[2]  for  example  shows  that  x2(t)  =  -0.125  Cos(t)  +  0  188  Cos(t)  -0.0521 
Cos(t). 

Cosine  Coefficients  through  order  25: 

c[0] 

{2.} 
c[2] 

{-0.125,  0.188,  -0.0521} 

c[4] 

{0.00594,  -0.0306,  0.0392,  -0.0194,  0.00298} 

c[6] 

{0.000976,  0.00137,  -0.00816,  0.0107,  -0.00652,  0.00186,  -0.000192} 

c[8] 

{-0.00022,  0.000518,  0.000446,  -0.00234,  0.00317,  -0.00221,  0.00085,  -0.000168, 

0.0000131} 
c[10] 
{-0.0000149,  -0.0000889,  0.000173,  0.000129,   -0.000701,  0.000984,  -0.000756, 

0.000351,  -0.0000976,   0.0000147,  -9.18xl0~7} 
c[12] 
{0.000011,  -0.0000129,  -0.000034,  0.000059,  0.0000364,  -0.000214,  0.000313, 

-0.00026,  0.000139,  -0.0000479,  0.0000104,  -1.26x  10~6,  6.56xl0"8} 
c[14] 

;-3.8 

0.000101,  -0.00009,  0.0000533,  -0.0000216,  5.91x10"°,  -1.04x10°,  1.06xl0"7, 
-4.7 
c[16] 


{-3.84xl07,  5.51x10°,  -4.1x10"°,  -0.0000121,  0.0000203,  9.83x10"",  -0.0000661, 
-4.75xl0"9} 
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{-5.56x10"',  1.1x10"',  2.13x10"°,  -1.4x10'°,  -4.24x10"°,  7.02*  10"°, 
2.46xl0"6,  -0.0000206,  0.000033,   -0.0000312,  0.0000201,  -9.22xl0"6,  3.xl0"6, 
-6.79xl0"7,  l.OlxlO"7,  -8.88x10  9,  3.47xl0"10} 

c[18] 

{8.17xl0"8,  -3.22xl0"7,  -8.11xl0-9,  7.85xl0"7,  -4.79x  10"7,  -1.46xl0"6, 
2.43xl0~6,  5.36xl0"7,  -6.42xl0"6,  0.0000108,  -0.0000108,  7.52xlO"6, 
-3.81xl0"b,  1.42xl0"6,  -3.85xl0"7,  7.42xl0"8,  -9.55x  10"9,  7.35x  10"10,  -2.55xl0"n} 

c[20] 

{2.46xl08,  2.84xl0"8,  -1.22xl0"7,  -1.12xl0"8,  2.85xl0-7,  -1.65x10  7,  -4.96xl0"7, 
8.45xl0"7,  8.07xl0"8,  -2.01xl0"6,  3.57xl0"6,  -3.76xl0"6,  2.78xl0"6,  -1.53xl0"6, 
6.37xl0"7,  -2.xl0"7,  4.65xl0"8,  -7.78xl0"9,  8. 82x10" 10,  -6.05xl0"n,  1.89X10"12} 

c[22] 

{-7.83xl09,  1.68xl0"8,  1.41xKT8,  -4.49xl0"8,  -6.79x  10"9,  1.02xl0"7,  -5.82xl0"8, 
-1.67x10  7,  2.94xl0^7,  -7.04xl0"9,  -6.26xl0"7,  1.18xl0"6,  -1.3xl0"6, 
1.02xl0"6,  -6.05xl0"7,  2.76xl0"7,  -9.72xl0"8,  2.63xl0"8,  -5.35xl0"9, 
7.91xl0"10,  -8.01xlO"U,  4.95xl0-12,  -1.41xl0"13} 

c[24] 

{-7.03xl0-10,  -3.68xl0"9,  6.03xl0-9,  5.84xl0"9,  -1.63xl0"8,  -3.2xl0"9, 
3.61xl0"8,  -2.08xl0~8,  -5.57xl0"8,  1.02xl0"7,  -1.41xl0"8,  -1.94xl0"7, 
3.89xl0"7,  -4.51xl0"7,  3.74xl0"7,  -2.36xl0"7,   1.16xl0"7,  -4.52xl0"8, 
1.38xl0"8,  -3.28xl0"9,  5.92xl0"10,  -7.84xl0"n,  7.17xlO-12,  -4.03xl0"13,  1.05xl0~14} 

Sine  Coefficients  through  order  25: 
s[l] 

{0.75,  -0.25} 

s[3] 

{-0.0273,  0.082,  -0.0608,  0.0122} 

s[5] 

{-0.00293,  -0.00879,  0.0199,  -0.0167,  0.00607,  -0.00075} 

s[7] 

{0.000743,  -0.000391,  -0.00254,  0.00555,  -0.00507,   0.0024,  -0.000562,  0.00005} 

s[9] 

{0.000039,  0.000262,  -0.000138,  -0.000759,  0.00164,  -0.00162,  0.000902, 

-0.000292,  0.0000499,  -3.46xl0"6} 
s[ll] 
{-0.000035,  -9.57xl0"6,  0.0000953,  -0.0000478,  -0.000227,  0.000505,  -0.000527, 

0.000332,  -0.000132,  0.000032,  -4.32x10"°,  2.45xl0"7} 
s[13] 
{1.51xl0~6,  -0.0000106,  -4.09xl0"6,  0.0000335,  -0.0000174,  -0.0000683,  0.000158, 

-0.000175,  0.000121,  -0.0000556,  0.000017,  -3.3xl0"6,  3.67xl0"7,  -1.76xl0"8} 
s[15] 
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{1.72xl(T6,  1.56  xlO-6,  -4.06x10°,  -1.73x10°,  0.0000115,  -6.5x10^, 

-0.0000206,  0.0000501,  -0.0000583,  0.0000437,  -0.0000225,   8.14xl0'6, 
-2.02xl0~6,  3.26xl0"7,  -3.08xKT8,  1.28xl0'9} 

s[17] 

{-2.66xl0"7,  4.13xl0~7,  6.46x10  7,  -1.46xl0"b,  -6.44xl0~7,  3.94xl0"6, 

-2.44x10 -6,  -6.16xl0~6,  0.000016,  -0.0000196,  0.0000157,  -8.91xl0'6, 

3.66xl0~6,  -1.09xKT6,  2.26xl0"7,  -3.12xl0"8,  2.56xl0~9,  -9.4xHrn} 

s[19] 

{-7.43xl0"8,  -1.34xl0"7,  1.57xHT7,  2.59xl0"7,  -5.16xl0"7,  -2.24xHT7,  1.35xl0~6, 
-9.17xl0'7,  -1.83xKT6,  5.15xl0"'5,  -6.62xl0"6,  5.62xlO"6,  -3.45xl0"6, 

1.58xl0"6,  -5.37xl0~7,  1.35xl0"7,  -2.41xl(T8,  2.91xl0'9,  -2.11xl0"10,  6.94xl0"12} 

s[21] 

{2.47xKT8,  -1.02xl0~8,  -5.54xl0"8,  5.45xl0"8,  9.78xl0~8,  -1.81xl0"7, 


-7.49xl0"8,  4.59xl0"7,  -3.43xl0^7,  -5.37xl0"7,  1.66x10°,  -2.24xl0"6, 

2.01xl0"6,  -1.32xl0"6,  6.59xl0"7,  -2.51xl0"7,  7.3xl0~8,  -1.59x  10"8, 
2.49xl0"9,  -2.66xl0~10,  1.73xl0-11,  -5.15xl0"13} 

s[23] 

{2.02xl0-9,  9.52xl0'9,  -3.45xl0"9,  -2.18xl0"8,   1.86xl0~8,  3.59xl0~8, 

-6.37xl0"8,  -2.4xl0~8,  1.56xl0"7,  -1.28  xlO"7,  -1.54xl0~7,  5.34xl0"7, 

-7.57xKT7,  7.15xl0"7,  -5.xl0~7,  2.69xl0"7,  -1.13xl0~7,  3.69xl0"8, 

-9.34xl0"9,  1.79xl0~9,  -2.5xl0"10,  2.4xl0-11,  -1.41x  10"12,  3.84x  10~14} 

s[25] 

{-1.86xl0"9f  -4.12xl0"10,  3.91xl0~9,  -8.87xl0"10,  -8.2xl0"9,  6.36xl0~9,  1.29x 
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-2.23xl0"8,  -7.37xl0~9,  5.29xl0'8,  -4.73xl0"8,  -4.31xl0~8,  1.72xl0"7, 

-2.56xl0"7,  2.54xl0"7,  -1.87xl0~7,  1.08xl0~7,  -4.9xl0~8,  1.77xl0"8, 

-5.08xl0"9,  1.14xl0"9,  -1.94xl0-10,  2.45xl0~n,  -2.14xl0-12,  1.15x  10"13,  -2.88xl0"15} 


APPENDIX  D  :  SOLUTIONS  TO  THE  VARIATIONAL  EQUATION 


This  appendix  includes  the  solutions  to  the  variational  equation  through  e  for  the  cases 
covered  in  Chapter  5  Note  that  these  solutions  can  be  used  to  check  if  the  code  in 
Appendix  C  is  working  correctly. 

Curves  out  of  Ao=0: 

(    Cos[t]       1  3Sin[t]       1  \ 

Lower  =  Cos[t]  +  Sin[t]  +  e +  -  Cospt]  + Sm[3t]  + 

v        8  8  8  8/ 

-/     19Cos[t]       1  5  7Sin[t]        1  5  \      w  5Cos[t] 

-2\ —  +  -Cospt] Cos[5t] +  — Sinpt] Sin[5t]  +  e3  — 

I         192  8  192  192         16  192  )        I 


e 
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19  7Cos[7t]       55Sin[t]        53  89Sm[5t]       7Sin[7t]  , 

Cospt]  + Cospt] — + Sinpt] + + 

256  768  1152  2304         768  2304  1152 

.1  4249Cos[t]       57Cos[3t]       3371Cos[5t]       2615Cos[7t]       61Cos[9t]       439Sin[t] 

e + +  +  - 


V     552960  2048  110592  221184  40960  552960 

107Sin[3t]       1631Sin[5t]       1885Sin[7t]       61Sin[9t]\      ,/1013Cos[t] 

+ +  +  e      + 


H 


9216  110592  221184  40960     ;        I    7372800 

6707  Cospt]       40205  Cos[5t]       180211  Cos[7t]       40709  Cos[9t]       5533  Cos[  lit] 
2211840  5308416  26542080  14745600  14745600 

8311Sin[t]       21127Sin[3t]       91811  Sin[5t]       319577  Sin[7t]       158179Sin[9t]       5533  Sin[  lit] 

+ 


22118400     2211840      5308416      26542080      44236800      14745600 
fCos[t]   1        3Sin[t]   1 


/Cost         1  JSmt         1  \ 

Upper  =  Cos[t]  -  Sin]t]  +  e -  Cospt]  + -  Sinpt]   + 

v      8  8  8  8/ 


,i     19Cos[t]       1  5  7Sin[t]        1  5  \      ,/     5Cos[t] 


ez +  -Cos[3t]- Cos[5t]+ — Sinpt]  + Sm[5t]  +  e  \- + 

I         192  8  192  192         16  192  j        I       576 

7                      19                    7Cos[7t]       55Sm[t]        53  89Sin[5t]       7Sin[7t] 
Cospt]- Cos[5t]  + + Sinpt]-      — — —  + 


256  768  1152  2304         768  2304  1152 

A  4249Cos[t]       57 Cospt]       3371  Cospt]       2615Cos[7t] 
6  I     552960      "       2048       +       110592       "       221184 

61Cos[9t]       439Sin[t]       107Sin[3t]       1631Sin[5t]       1885Sin[7t]       61Sin[9t]\ 
40960     ~     552960  9216  110592      +       221184  40960    ) 

5{     1013Cos[t]       6707Cospt]       40205Cos[5t]       180211  Cos[7t]       40709Cos[9t]       5533Cosfllt] 


I       7372800            2211840              5308416               26542080  14745600  14745600 

8311Sin[t]       21127Sin[3t]       91811Sin[5t]       319577Sin[7t]  158179Sin[9t]  5533 Sin[  1 1 1]  \ 

22118400  "       2211840  +      5308416  ~       26542080      +  44236800  14745600    j 

Curves  out  of  A0=3/2: 


63 

,/     71  23  7  \ 

Lower  =  Cos[2t]  +  e2\ Cos[2t]  + Cos[4t] Cos[6t]    + 

1     576  144  192  J 

,i  15443Cos[2t]       2809Cos[4t]       1427Cos[6t]       1127Cos[8t]       119Cos[10t] 

-  + + 


1658880  82944  36864  69120  552% 

1  \      %{    73Sin[2t]        47     .  239Sm[6t]         5 


—  Sin[2t]-  -Sin[4t]   +e3 + Sin[4t] + Sin[8t]  + 

1 24  6  J        I        1728  576  4608  576  ) 

5(    5407Sin[2t]        1 12483  Sin[4t]       285299Sin[6t]       65987  Sin[8t]       22187  Sin[10t]       1007Sm[12t1  \ 
i"     13271040  9953280       +       13271040  4147200      +       4423680  1843200     J 


11  1 

Upper  =  e  |  -  —  Cos[2t]  +  —  Cos[4t]  |  + 


3(55Cos[2t]       181Cosf4t]       239Cos[6t]       5Cos[8t]  \      5< 
e  I      1536  3456  9216  1152     )  +  €  \ 


319213  Cos[2t] 
79626240 


218417Cos[4t]       378463Cos[6t]       1117Cos[8t]       22187Cos[10t]       1007Cos[12t]  \ 


+ 


19906560      26542080      129600       8847360       3686400 


V 


23  7  \ 

Sin[4t] Sin[6t] 

288  384  J 


1  ,/     121Sin[2t] 
—  Sm(2t]  +  c + 

2  I         1152 

4(  46007 Sin[2t]       4385Sin[4t]       4883Sin[6t]       1127Sin[8t]       119Sin[10t] 
3317760      "       165888      +       221184      "       138240      +       110592 


Curves  out  of  A0=4: 

Lower  =  Cos[3t]  -  Sin[3t]  + 

/3Cos[t]        3  3  3Sin[t]       359  3  \      it    329Cos[t] 

e — Cos[3t]-  —  Cos[5t]  + + Sin[3t]  -  —  Sin[5t]   +  e 

I       8  16  16  8  240  16  7        I         640 


+ 


233Cos[3t]       479Cos[5t]        27  51Sin[t]       73913Sin[3t]        15  27 

+ Cos[7t]- Sin[5t]  + Sin[7t][  + 

1280  1280  640  128  57600  256  640 

w89813Cos[t]       76151Cos[3t]       139  67Cos[7t] 

e Cos[5t]- + 

I      153600  230400  600  2048 

469Cos[9t]       831Sin[t]       82393Sm[3t]       679Sin[5t]       5311Sm[7t]       469Sm[9t]  \ 

46080       +      2048  57600  5120  51200  46080     )  + 

4(    385217  Cos[t]       105987247  Cos[3t] 

E  I"       614400       +         309657600 

389277  Cos[5t]       30401Cos[7t]       327931  Cos[9t]       547Cos[llt]       273059  Sin[t] 

1638400  409600  11059200  215040  614400 

37616912581  Sin[3t]       18809Sin[5t]       249571  Sm[7t]       27691  Sin[9t]       547Sin[llt]\ 

23224320000  327680  6144000       +       2211840  215040     )+ 

5/  4244707848 lCos[t]       33657149027  Cos[3t]       9222921259 Cos[5t]       1330129 Cos[7t] 

I       6193152(XKX)  92897280000  30965760000  39321600 

46069577  Cos[9t]       88301  Cos[l  It]       89371  Cos[13t]       400303339  Sin[t]       82024682603  Sn[3t] 

9289728000       +       20643840  137625600  825753600  46448640000 

13130773 Sin[5t]       160571  HSin[7t]       1 1346613 Sin[9t]       883481  Sm[ lit]       89371  Sin[  13 1]  \ 


+ 


2064384a)      589824000      371589120      103219200      137625600 


■ 


64 


Upper  =  Cos[3t]  +  Sm[3t]  + 

I     3Cos[t]        3  3  3Sin[t]       359  3  \      ./     329Cos[t] 

e +  — Cos[3t]+  — Cos[5t]  + +  Sin[3t] Sin[5t]    +e2 

I         8  16  16  8  240  16  V        V         640 


233Cos[3t]       479Cos[5t]        27  51Sin[t]       73913  Sin[3t]        15  27 

■  + 


1280  1280  640 

89813Cos[t]       76151  Cos[3t]       139  67Cos[7t] 


51Sm|t|        VjyiJSmlJtJ         15  LI  \ 

Cos[7t]  + + + Sin[5t] Sin[7t] 

128  57600  256  640  ) 


i     HVSUCoslt]        76151  Cos|3t]        139 

+ + Cos[5t]  + 

V    153600      230400     600 


2048 

469Cos[9t]       831Sm[t]       82393  Sin[3t]       679Sin[5t]       5311Sin[7t]       469Sin[9t] 

■  + + + + 


46080  2048  57600  5120  51200  46080 

3852 17  Cos[t]       105987247  Cos[3 1] 

614400       +         309657600         + 
389277  Cos[5t]       30401  Cos[7t]       327931  Cos[9t]       547Cos[llt]       273059Sin[t] 


1- 


1638400  409600  11059200  215040  614400 

37616912581  Sin[3t]       18809Sin[5t]       249571  Sin[7t]       27691  Sin[9t]       547Sm[llt] 

■  + 


23224320000       327680      6144000      2211840      215040 


42447078481  Cos[t]   33657149027  Cos[3t]   9222921259Cos[5t]   1330129  Cos[7t] 
e + + + 


61931520000       92897280000       30965760000       39321600 
46069577  Cos[9t]   88301  Cos[llt]   89371Cos[13t]   400303339  Sin[t]   82024682603  Sin[  3 1] 


+ + + + 


9289728000      20643840      137625600      825753600       46448640000 
13130773  Sin[5t]   160571 1 1  Sin[7t]   1 1346613  Sin[9t]   883481  Sin[  lit]   89371  Sin[13t] 


2064384O)      589824000      371589120      103219200      137625600 


APPENDIX  E  :  STABILITY  TRANSITION  CURVES 


This  appendix  includes  the  stability  transition  curves  found  in  Chapter  5.  Curves  are  given 
analytically  through  order  10  and  numerically  through  order  30. 

Curves  out  of  A0=0: 


c   c2   c3   7c4   3e5 

247  c6 

11657e7 

1224811c8 

2   8   32   384   1024 
65987129c9    248740367c10 

442368 

21233664 

2548039680 

611529523200   36691771392000 

Lower  eA  = 

-05e  +  0125c2  +  0.03125e3  -  0.018229c4  +  0.0029297e5  - 
0.00055836e6  -  0.00054899c7  +  0.00048069c8  -  0.00010791  e9  +  6.7792x  10-6  e10  + 
0.000023325  c11  -0.000020161c12  +3.8825  x  KT6e13  +  7. 1301  x  KT7e14-  1.1828x  10_6e15  + 
8.9935 x  10_7e16  -  1.1027x  KT7e17  -  99766x  10_8e18  +  6.0926x  KT8e19  -  3.6467  x  KT8e20  + 
1.5256x  KT11  e21  +  8.7439x  KT9e22  -  2.8529x  10"9  e23  +  98924x  KT^e24  +  3.551  x  10~10  e25  - 
6.2361  x  10"10e26  +  9.9193 x  10_11e27  +2.5033 x  10"11  e28-  3.8338x  10" 1 1 e29  +  3.6974 x  KT11  e30 


e  e2       e3       7e4   3e5 

247  e6    11657  c7 

1224811c8 

OCTtA  _    +     —     —      — 

v              2   8   32   384   1024 

442368  +  21233664  4 

2548039680 

65987129c9    248740367c10 

611529523200  +  36691771392000 

Upper  eA  = 

05e  +  0.125e2  -  0.03125c3  -  0.018229c4  -  0.0029297c5  - 
0.00055836c6  +  000054899c7  +  0.00048069c8  +  0.00010791  e9  +  6.7792x  10"6  e10  - 
0.000023325  c11  -0.000020161c12  -  3.8825  x  10~6e13  +  7.1301  x  10-7e14  +  1.1828x  10"6e15  + 
8.9935x  KT7e16  +  1.1027x  10-?e17  -  99766x  10"8  e18  -  6.0926x  10~8e19  -  3.6467x  lO"8  e20  - 
1.5256x  10"11  e21  +  8.7439x  10"9e22  +  2.8529x  10"9 e23  +  98924x  lO-^e24  -  3.551  x  10"10 e25  - 
6.2361  x  lO-^e26-  9.9193  x  lO^e27  +2.5033  x  10" 1 1  c28  +  3.8338 x  10"nc29  +  3.6974 x  10"11  c30 

Curves  out  of  A0=3/2: 

109c4        3419c6        297305c8         7797193367c10 


LowereA  = + 


2       6        3456        1990656      573308928       82556485632000 
Lower  eA  = 


66 
1.5-  0.16667 e2  +  0031539c4  -  0.0017175e6  -  O.0OO51858e8  +  0.(X)0094447 f 10  + 
0.000013 172  e12-  5.7415 x  10"6e14-  1.1848x  10"7e16  +  3.354 x  10"7e18  -  2.8733 x  Hr8e20- 
1.7574  x  UT8e22  +  3.7778  x  KTV4  +7.2633 x  10" '"f26-  3.3523 x  K)-10e28-  1.1097x  10"11  e30 

3       e2       53 e4        6983 e6  740213 e8  6745666577 e 10 

Upper  e  A  =  —  +  —  - 


2       3       3456       1990656      2866544640       82556485632000 


Upper  eA  = 

1.5+  0.33333  e2-  0.015336e4-  0.0035079e6  +  0.00025822 e8  +  0.00008171  e10  - 
0.00001 5556 e12-  1.5903 x  10~6e14+  1.0185x  lO'V6  -  3.4918 x  Kr8e18  -  5.8479 x  K)-8e20  + 
89289x  lO^e22  +  2.7575x  lO^e24  -  9.066x  lO^e26  -  8.0056 x  10-ne28  +  7. 1452 x  10"11  e30 

Curves  out  of  A0=4: 

9e2       5e3       2877 e4       2537 e5        54707 e6         1 1663959 e7 
Lower  eA  =4 + + + 


32        64        40960        737280       15728640       14155776000 
17260922723 e8  43486217423 e9  2804162996941 e10 


12683575296000       1223059046400000       1461 1478740992000 

Lower  eA  = 

4  -  0.28125e2  -  0078125e3  +  O070239e4  +  0.003441  e5  -  0.0034782e6  + 
0.00082397e7-0.O013609€8  -0.000035555 e9  +  0.000 19 1 92 e10-  0.0000161 18c11  + 
0.000037433  e12  +  1.4085  x  KT6e13  -  0.00001 1912 e14  +  2.053  x  10_7e15  -  6.2434  x  KT7e16  - 
6.81 13 x  10"8e17  +  7.0367 x  10"7e18  +  1.5598x  10"8e19  -  4.6186 x  10~8 e20  +  2.0464 x  10"9e21  - 
3.6947  x  10"8e22-  2. 1736  x  lO^e23  +  7.4707  x  \0~9  e24  +  9323  x  10"  Ue25  +  1.5107  x  10"9e26  + 
1.8088x  KT10e27  -  6.9398x  KT10e28  -  2.4687x  KT11  e29  -  2.0128x  KT11  e30 

9e2       5e3       2877e4       2537e5        54707e6         11663959e7 


UppereA  =4 + + 


32    64    40960   737280   15728640   14155776000 
17260922723 e8    43486217423 e9    2804162996941 e10 


12683575296000   1223059046400000   1461 1478740992000 

Upper  eA  = 

4  -  028125e2  +  0078125e3  +  0070239e4  -  0.003441  e5  - 0.0034782e6  - 
0.00082397  e7  -  0.0013609e8  +  0.000035555  e9  +  0.00019192e10  +  0.000016118en  + 
0.000037433 e12  -  1.4085 x  KT6e13  -  0.00001 1912 e14  -  2.053 x  HT7e15  -  6.2434x  10-7e16  + 
6.8113x  10"8e17+  7.0367  x  KrV8  -  1.5598x  10"8f19-  4.6186 x  10-8e20-  2.0464 x  10"9e21  - 
3.6947  x  10-8e22  +  2. 1736 x  lO^e23  +  7.4707 x  10"9e24-  9.323  x  10" ne25  +  1.5107 x  HT9  e26  - 
1.8088x  10-10e27  -  6.9398x  K)-10e28  +  2.4687x  KT11  e29  -  2.0128x  KT11  e30 


on     51NPS  i !  TOM 
,/gg  22527-15"? 


